{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "505cad3a-c73f-4500-97e6-67afc97fd775",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>p</th>\n",
       "      <th>key</th>\n",
       "      <th>burglary</th>\n",
       "      <th>earthquake</th>\n",
       "      <th>alarm</th>\n",
       "      <th>john</th>\n",
       "      <th>mary</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.950</td>\n",
       "      <td>alarm</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.940</td>\n",
       "      <td>alarm</td>\n",
       "      <td>1</td>\n",
       "      <td>-1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.290</td>\n",
       "      <td>alarm</td>\n",
       "      <td>-1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.001</td>\n",
       "      <td>alarm</td>\n",
       "      <td>-1</td>\n",
       "      <td>-1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.900</td>\n",
       "      <td>john</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.050</td>\n",
       "      <td>john</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.700</td>\n",
       "      <td>mary</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.010</td>\n",
       "      <td>mary</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>-1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.001</td>\n",
       "      <td>burglary</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.002</td>\n",
       "      <td>earthquake</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       p         key  burglary  earthquake  alarm  john  mary\n",
       "0  0.950       alarm         1           1      0     0     0\n",
       "1  0.940       alarm         1          -1      0     0     0\n",
       "2  0.290       alarm        -1           1      0     0     0\n",
       "3  0.001       alarm        -1          -1      0     0     0\n",
       "4  0.900        john         0           0      1     0     0\n",
       "5  0.050        john         0           0     -1     0     0\n",
       "6  0.700        mary         0           0      1     0     0\n",
       "7  0.010        mary         0           0     -1     0     0\n",
       "8  0.001    burglary         0           0      0     0     0\n",
       "9  0.002  earthquake         0           0      0     0     0"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# ===== 代码块1：导入库和加载数据 =====\n",
    "import os\n",
    "import random\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import networkx as nx\n",
    "\n",
    "burglary_pd = pd.read_csv(\"./burglary_cpts.csv\")\n",
    "burglary_pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "bfcf3a24-5cc3-4704-853e-135ef9b6edc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ===== 代码块2：定义基础类 =====\n",
    "import os\n",
    "import random\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import networkx as nx\n",
    "\n",
    "class Thing:\n",
    "    def __init__(self, name=''):\n",
    "        self.name = name\n",
    "    \n",
    "    def __repr__(self):\n",
    "        return self.name\n",
    "\n",
    "# 贝叶斯网络的节点结构是向上的\n",
    "class Node(Thing):\n",
    "    def __init__(self, name, parents=None, loc=None, state=-1):\n",
    "        super().__init__(name)\n",
    "        self.parents = parents if parents is not None else {}\n",
    "        self.loc = loc\n",
    "        self.state = state\n",
    "    \n",
    "    def add_parent(self, parent_name, parent_weight):\n",
    "        self.parents[parent_name] = parent_weight\n",
    "\n",
    "class Graph(Thing):\n",
    "    def __init__(self, name):\n",
    "        super().__init__(name)\n",
    "        self.node_list = []\n",
    "        self.node_map = {}\n",
    "    \n",
    "    def add_node(self, node):\n",
    "        self.node_list.append(node)\n",
    "        self.node_map[node.name] = node\n",
    "    \n",
    "    def draw(self, font_size=10, k_layout=1, node_size=1000, title=''):\n",
    "        plt.rcdefaults()\n",
    "        plt.rcParams['figure.dpi'] = 150\n",
    "        g = nx.DiGraph()\n",
    "        edges = []\n",
    "        edges_map = {}\n",
    "        \n",
    "        for node in self.node_list:\n",
    "            for key, val in node.parents.items():\n",
    "                edges.append((node.name, key, 1/val))\n",
    "                edges_map[(node.name, key)] = val\n",
    "        \n",
    "        g.add_weighted_edges_from(edges)\n",
    "        \n",
    "        pos_map = {}\n",
    "        nodes_label = {}\n",
    "        for node in self.node_list:\n",
    "            if node.loc is not None:\n",
    "                pos_map[node.name] = node.loc\n",
    "            nodes_label[node.name] = node.name\n",
    "        \n",
    "        rand_seed = random.randint(1, 1000)\n",
    "        if len(pos_map) > 0:\n",
    "            locs = nx.spring_layout(g, pos=pos_map, fixed=pos_map.keys(), seed=rand_seed, k=k_layout)\n",
    "        else:\n",
    "            locs = nx.spring_layout(g, seed=rand_seed, k=k_layout)\n",
    "        \n",
    "        color_arr = []\n",
    "        for node in self.node_list:\n",
    "            color_arr.append(node.state)\n",
    "        \n",
    "        nx.draw_networkx_nodes(g, pos=locs, node_size=node_size, node_color=color_arr, cmap=plt.cm.cool)\n",
    "        nx.draw_networkx_labels(g, pos=locs, labels=nodes_label, font_size=font_size, font_color='black')\n",
    "        nx.draw_networkx_edges(g, pos=locs, node_size=node_size)\n",
    "        plt.title(title)\n",
    "        plt.show()\n",
    "# ===== 代码块3：定义贝叶斯网络相关类 =====\n",
    "class BayesNode(Node):\n",
    "    def __init__(self, name):\n",
    "        super().__init__(name)\n",
    "        self.cpt = {}\n",
    "    \n",
    "    def cpt_score(self, evidence):  # 只返回为1的概率\n",
    "        # 我们假设，调用这个方法的时候，该节点的所有支撑变量都进了evidence(这需要推理过程确保)\n",
    "        key_tuple = []\n",
    "        for parent in self.parents:\n",
    "            if parent in evidence:\n",
    "                key_tuple.append(evidence[parent])\n",
    "        key_tuple = tuple(key_tuple)\n",
    "        # 对于先验概率，其tuple应该是一个空值，这样可以直接算出需要的p\n",
    "        p = self.cpt[key_tuple]\n",
    "        return p\n",
    "    \n",
    "    def __repr__(self):\n",
    "        return \"{}{}\".format(self.name, tuple(p for p in self.parents))\n",
    "\n",
    "class BayesNet(Graph):\n",
    "    def __init__(self, name='', cpt_pd=\"data/burglary_cpts.csv\"):\n",
    "        super().__init__(name)\n",
    "        self.cpt_pd = pd.read_csv(cpt_pd)\n",
    "        self.parse_cpts(self.cpt_pd)\n",
    "    \n",
    "    def parse_cpts(self, data_pd):\n",
    "        # 首先加载各个点\n",
    "        # 注意这里必须严格确保源文件中列的顺序，因为检查的时候是从先验概率一层一层补充检查的\n",
    "        for col in data_pd.columns[2:]:\n",
    "            self.add_node(BayesNode(col))\n",
    "        \n",
    "        # 构造各个点内部关系\n",
    "        for idx, row in data_pd.iterrows():\n",
    "            p, key = row['p'], row['key']\n",
    "            row['p'], row['key'] = 0, 0\n",
    "            row = row.where(row != 0).dropna()\n",
    "            row_node = self.node_map[key]\n",
    "            for col, val in row.items():\n",
    "                row_node.parents[col] = 1\n",
    "            vals = tuple(val for _, val in row.items())\n",
    "            row_node.cpt[vals] = p\n",
    "    \n",
    "    def ask(self, x, evidence, print_formula=False):\n",
    "        result = {1: 0, -1: 0}\n",
    "        if x in evidence:  # 如果x的状态已知\n",
    "            result[evidence[x]] = 1.0\n",
    "        else:  # 拿出来所有的node，逐一排查\n",
    "            cal_seq = []\n",
    "            cal_seq.append(\"1:\")\n",
    "            result[1] = self.recurse_ask(self.node_list, {**evidence, x: 1}, cal_seq)\n",
    "            cal_seq.append(\"\\n-1:\")\n",
    "            result[-1] = self.recurse_ask(self.node_list, {**evidence, x: -1}, cal_seq)\n",
    "            result_sum = result[1] + result[-1]\n",
    "            result[1] = result[1] / result_sum\n",
    "            result[-1] = result[-1] / result_sum\n",
    "        \n",
    "        if print_formula:\n",
    "            print(\"\".join(cal_seq))\n",
    "        return result\n",
    "    \n",
    "    def recurse_ask(self, check_list, evidence, cal_seq):\n",
    "        if not check_list:\n",
    "            cal_seq.append(\"1\")\n",
    "            return 1.0\n",
    "        \n",
    "        check_node, rest = check_list[0], check_list[1:]\n",
    "        if check_node.name in evidence:  # 如果这个状态是一个已知的值\n",
    "            cpt_score = check_node.cpt_score(evidence)\n",
    "            if evidence[check_node.name] == -1:\n",
    "                cpt_score = 1 - cpt_score\n",
    "            cal_seq.append(\"({}*\".format(cpt_score))\n",
    "            recurse_val = self.recurse_ask(rest, evidence, cal_seq)\n",
    "            cal_seq.append(\")\")\n",
    "            p = cpt_score * recurse_val\n",
    "        else:  # 如果这个状态未知，就需要计算+-1的所有情况\n",
    "            cpt_score = check_node.cpt_score(evidence)\n",
    "            cal_seq.append(\"({}*(\".format(cpt_score))\n",
    "            p = cpt_score * self.recurse_ask(rest, {**evidence, check_node.name: 1}, cal_seq)\n",
    "            cal_seq.append(\")+{}*\".format(1 - cpt_score))\n",
    "            p += (1 - cpt_score) * self.recurse_ask(rest, {**evidence, check_node.name: -1}, cal_seq)\n",
    "            cal_seq.append(\")\")\n",
    "        return p\n",
    "\n",
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6e5dd8ad-4e11-4e10-8958-76a368006918",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwYAAAJICAYAAADb4/N8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAABcSAAAXEgFnn9JSAABwE0lEQVR4nO3dCZyNdfvH8e9g7DtZipCoqOy7LNmzpMVa0qJFSkqrkh75p6RFepRQ6qFQ2ffsIftSKdkiRSn7Ppb5v67fOTMZZhhjZu6zfN6v13lmzrnnnPObo4f7un/XEhEdHR0tAAAAAGEtjdcLAAAAAOA9AgMAAAAABAYAAAAACAwAAAAAEBgAAAAAMAQGAAAAAAgMAAAAABAYAAAAACAwAAAAAGAIDAAAAAAQGAAAAAAgMAAAAABAYAAAAADApEvqx1CgQAEdPnxYV155JZ8kAAAA4LHffvtNWbJk0Z9//pm6OwYWFJw4cSKpTwcAAACQjOzc3M7RU33HIGanYN26dUl+cwAAAADJo3Tp0pf0fGoMAAAAABAYAAAAACAwAAAAAEBgAAAAAMAQGAAAAAAgMAAAAABAYAAAAACAwAAAAACAITAAAAAAQGAAAAAAgMAAAAAAgKR0Xi8gmByW9Kuko/7bcUnpJWWSlFFSEUk5vF4kAAAAkAQEBucJAtZIWnnG7WdJpy/wvKslVZRUwX8rT7AAAACAIEBgcIYjkkZLGixpeSKCgPhs8t9GnfHY9ZI6SeooKWcyrhcAAABILtQYSNooqbukQpLul7Q0iUFBQn6U1E3SFZIe8u9EAAAAAIEkrAODmZIaSSop6W1Je1NhR2KIpHKSqvt3J6JT+D0BAACAxAjLwGCPpA7+oMCCAy98J6mtpFsk/e7RGgAAAICwDQwmSSotaYQCw3T/ej5m9wAAAAAeShNuuwQtJP2pwHJA0gPsHgAAAMBDYREYLPV3BgqUXYIL7R5M8XohAAAACDshHxjMllRP0k4FB9s9aCnpC68XAgAAgLAS0oHBRH96jg0rCyYnJd3ln6cAAAAApIaQDQy+kdRKUpSCkxUiPyLpM68XAgAAgLAQkoHBd/50nGANCs5kA9fGe70IAAAAhLyQCwys41Bz/zCxUHDKP+9gndcLAQAAQEgLqcDA0m86S9qt0HJc0r3+2gMAAAAgJYRUYDAqhNNuVkjq7/UiAAAAELLShFIK0WMKbb1IKQIAAEAKSRNKKUQ23TiUWTE1KUUAAABICSERGHwZwilE8aUUveP1IgAAABBy0oTCbkFfhZe3QqQVKwAAAAJH0AcGSyStUXj5S9I4rxcBAACAkBL0gcEghadw/b0BAACQMoI6MPhb0hiFpwWSfvB6EQAAAAgZQR0YfBzmufYfeL0AAAAAhIygDQxOS/pQ4e1/kg56vQgAAACEhKANDDZJ2nqpL1KnjhQRIW295FfyxCFJS8+4/8orrygiIkLDhw/3cFUAAAAIRkEbGKxMzA8VLeo78Q9hifocAAAAgLAODMJk4BkAAAAQtoEBJ8Q+BEgAAAAI+MBg+/bteuyxx1S8eHFlzJhRuXPnVrNmzbR48eI4PxcdHa0vvvhCbdu2VcmSJZUlSxZly5ZNlStX1qBBg3T6tJUa/8vufffKK740IcunX7ZMatZMypPH99i77/q+btvme4J9H3Oz9KL4jB8vVa0qZcki5c4ttWsn/f57/D+7Z4/02GPS5ZdLGTNKpUpJAwbYLxL/e5y51sSmPNlrffGF1LatVLKkb13ZskmVK0uDBkn+z+RXW855/gxiPt9u3bq5+oOaNWtq3759scd+/vln3XvvvSpcuLAyZMig/Pnzuz+HdevWXeBVAQAAEErSpdQLf/fdd2ratKn27t2ra665xn3/999/a8aMGZo+fbpGjhypNm3auJ89fvy42rdvrzx58qhUqVIqX768du/e7QKILl26aNmyZXEKajed2aZ0wQLpoYd8J88NG0o7dkg33yx17Ch99ZV0+LDv+xh58567WDvRfvtt6aabpFtukZYulUaNklaulNaulTJl+vdn9+6Vata0M2qpQAHp1lt9gcLTT0ubbGXJ5PhxqX17X7BjgUf58tLu3ZIFVV26+IIh/2eySlL9BF7m5MmTuv/++/W///1PTZo00VdffaXMmTO7Y+PHj3dBgH3+ZcuWVdWqVV0wN2bMGE2aNEnTpk1TrVq1ku93AgAAQOCKTqJSpUq5W3z2798fXbBgwei0adNGjxgxIs6x5cuXR+fKlSs6a9as0bt27XKPnThxInrcuHHRUVFRcX7WjlesWDHaljl//vzYx8fbJfBevdzj7vbGG9HusvjZtyJFfMfjO2a32rV9xzNnjtbixf8+fvhwtKpX9x0bNizucx55xPd448a+n4t5fOnSaGXN6jtm73vmc2LW+skniV/niRPRGjcuWlFRcR/ftSta/s9E8+e7xwb4P5de/vf55JNP3P0jR45EN2vWzD3Wrl27OJ/vr7/+Gp0lSxb35/DNN9/E+dynTZsWHRkZGV24cOHo48ePJ/U/EQAAAKSi852fJ0aKpBJ9/PHH2rlzp0tfueuuu+Icq1ixonr27KlDhw5pxIgR7rF06dKpZcuWioyMjPOzl112mfr27eu+nzBhQuzjh8/8oRtukJ555tIW/OSTUrVq/963K+pPPfXvjkTsGx+WPv1USpNGev9938/FsBQfu5KfXNKlk1q2lM76THTZZZL/M5H/MzkSz9P379+vRo0aafLkyXr00UfdDs2Zn++7776rw4cPu8+3fv24+w2NGzdW586d3e7BlClTku93AgAAQHilEs2cOdN9vf322+M9fpOl7MiyYZbFeXzNmjXuudu2bdORI0fcBfKDB30jvDZu3Bj7c8fOfJLVFlxqS1JLQTqbpSaZnTv/fcxSi44e9QUBxYuf+xyrS3jjDSWrNWvsA/XVSxw54qs98H8m8n8mcT4PSbt27VKdOnXc5/nSSy/p1VdfTdKf0Xvvvef+jG677bbk/Z0AAAAQHoHBVv/AsBo1apz35/755x/3NSoqyhXAWgFyQmIChHNOhK+88lKXKxUqdO5jVugbk+sfw+oXTJEi8b9OQoXNSREVJd17r68AOSH+z+ToWQ+/+OKLrrbArvrHFxSc+Wd0xRVXJOrPCAAAAKEtRQKDmC5Cd955p+swlJBrr73WfX377bddUHDDDTeoX79+rvg4V65cLvVlw4YNrnjZdg/iXbR1BbpUlhrkpbO6LjlWDG1BgaVK9evnKz7OlcuXWrRhg3TNNb7dA0lnJRu5XQArLLaCYyvqtk5E576l7z07nlmYHY8qVapcym8GAACAcA4MChUqpF9++UXPP/+8KlSocMGfHzdunPtqwUHp0qXjHNuyZcs5P58MoUDSFCzo+xrTBvVsCT2ePr3v66FD5x47dUr6889zH/d/Ji44OOsz0Vmfydmfh3UfstoOC8zse+sEVb169XP+jDZv3qy33nrLdYMCAABAeEuRS+UNGjSIc8J/IdbSNOZk9WzWOvNsZzQPPb+YE/KTJ5UsLMix1qVWaxBPwOJanJ4voLAr/WebO1c6ceLcx/2fSbxpTmd9JvF9Hi1atHCfnbUitWLiJUuWXNKfEQAAAEJbigQGDz/8sPLly+fSgj766KNzBpRZ/rtdxf7xxx/dfRtqZj788MM4P2c99z/77LNzXv/yxC7EBpCZX35RssiaVerQwXeV//HHfYXIMVas8HUqik/MLADrwuTP7Xd+/VXq2jX+58QUP5/1mbjZDGd9Jv6w4xzW6WnUqFE6evSoCw7OLPbu3r27MmXKpKefflpjx44957kWUNjn/3tCQ94AAAAQUlIkMMiZM6drL5ojRw4XJBQtWlS33HKLS2+pV6+ea0NqJ6qb/APBnn32WaVNm9alHlk7U8uLr1Spklq1aqUnrZXoWcokdiEtWvi+1qvn6xjUqZP0/POX9stZq1DL75861deZyCYTN2rka3d6zz3xP8d+zo7ZLkDZsr51WYtQqx+4/vr4i5mffVZKm9a33ooVfcPOKlWSWrXytVc9Q8XzLNfqDSxFy1qTWvvSFRbASLr66qvd4ydOnNAdd9yhEiVKuF2Gdu3auaFmll5knz/FxwAAAOEhxapubYruDz/84E76s2fPrvnz57uCWGtFWrt2bTfJOKZ/vp2ILly4UDfffLOrKbDe++nTp9fXX3/tJh+fLasNME7MIuxq/Esv+a70f/21NGxYwuk+iZU7t7RokdS5s6/4d/x46bffpNdflwYOTPh5Q4b4TvKzZ5dmzPDtHLzwQsJdh2yXYeFC3xRnS1uaPNmXGmW/xxmfifVOuvoCS7ZaA5tjYJ2dGjZsqFWrbFayDW2+Vd9//72bcxAREaFvvvnGzS2wdqfNmzd3qUg2iRoAAAChL8KmnCXliTFFwuvWrZMX7Nr8/xSAbKaC7QCcmTKUgmpLmpcq7wQAAIBAdqnn5x736Uy6C/c6Cg98DgAAAEgOBAZBjs8BAAAAYR0YlLNupF4vIgAwfgwAAABhHRjYPOXWCkBWspFK9QU3WcOjVHknAAAAhLqgDQzMowpv4f77AwAAIPkEdWBQVVJZhad8NqPA60UAAAAgZAR1YBAhqbPC04PUWAAAACAZBXVgYNpLyq7w+0N7yOtFAAAAIKQEfWBgU5DvV3i5VdKVXi8CAAAAISXoAwPzsqTLFR6sG9NbXi8CAAAAISckAoNckj5SeOgnqZjXiwAAAEDICYnAwDSVdI9CWx1Jj3i9CAAAAISkkAkMzLshnFJkKUSPrlihlcuX6+TJk14vBwAAACEmnUJITEpRM4WeJ3fsUOtKldz3kZGRuuGGG9ztxhtvdDf7Pn/+/F4vEwAAAEEqpAKDmJSiJyW9o9Bxh6TuWbPq9XTp3G7BiRMntGrVKnc7U758+VyQ0KpVKz30EA1NAQAAEKapRDH6S+qo0NBQ0khJObNn1yuvvHLen921a5dmzZqlJ598UqdOnUq1NQIAACD4pQnVX2qopNsV3GpIGispg//+E088obx5817wed27d1fatGlTfH0AAAAIHSEZGMTkSI2SdJeCU31J0/1FxzGyZs2qZ5999rzPq127tnr16pXi6wMAAEBoCdnAwERK+kxSFwWX2yRN9k91Ptujjz7qagniky5dOr3++uvsFgAAAOCihXRgEPMLDpQ0TFJ2Bbb0kl6X9OUZ6UNny5IlS4K7BlaY3LhxY40dawlIAAAAQOKFfGBgIiTdL+lHSY0UmKwR6WpJz0m60PX+zp07n9Oa1HYSqlWrpv379+uOO+7Q448/ruPHj6fomgEAABA6wiIwiFFY0jR/YXL2ANslWCypVCKfkzlzZj33nIUQPtWrV9e7776r+fPnx+4mvP/+++7xTZs2pdDKAQAAEErCKjCI2T14wL970MLjtdQ6Y5fgYgdKPPLIIypRooQKFiyoUaNGuaFndnvjjTc0ZcoU5cmTx805KF++vEaPHp1CvwEAAABCRdgFBmfuHkyQtEHSUzYnIJXeN7OkByXZaLL5F7FLcLZMmTLp+++/dzsChQvbb/OvW265RWvWrNFNN92kgwcPqm3bti6QOHr0aLL8DgAAAAg9YRsYxCgh6S1Jf/gLlMun0PtcI2mA/30+klQuGV4zY8aMLq0oPoUKFdKcOXP04osvKiIiQoMHD1bVqlX1yy+/JMM7AwAAINSEfWAQI7O/QHmF/2bBQnv/Cb2lH12sYpLulNRX0gJJP0vqmoo7EzHtS/v06aMZM2a4Fqe2w1ChQgWNGDEiFVcBAACAYBARHR0dnZQnli5d2n1dt26dQt0BSWv8AcNGScckWVLOcX/xcCa7ei+pqKQK/l2HPAosO3fu1F133aW5c+e6+/fff78GDhyY4I4DAAAAgsulnp8TGISRU6dOuR2E//znP7I/9lKlSunLL790XwEAABDcLvX8nFSiMGITkXv16qXZs2erQIEC+umnn1SxYkV98sknLlAAAABA+CIwCEN169Z1XYsaNGjgOhVZWlHHjh116NAhr5cGAAAAjxAYhCmbnDx9+nT93//9n9KkSaP//e9/qlSpkitQBgAAQPghMAhjFhD06NFD8+bN0xVXXKH169erSpUq+uijj0gtAgAACDMEBnCD0Cy1qEmTJjp27JgefvhhtW/fXgcOWD8mAAAAhAMCAzh58+bV5MmT1a9fP1ekPGrUKDfzYPXq1V4vDQAAAKmAwABxUoueeeYZffvtt7ryyiu1adMmNy150KBBpBYBAACEOAIDnKNatWpup6BFixaKiopSly5d1Lp1a+3fv9/rpQEAACCFEBggXrlz59b48eP1zjvvKDIyUl999ZXKlSun5cuXe700AAAApAACAyQoIiJC3bp106JFi1S0aFH9+uuvqlGjht59911SiwAAAEIMgQEuyOYbWGrR7bffrhMnTujJJ5/Ubbfdpj179ni9NAAAACQTAgMkSs6cOV060fvvv6/06dNrwoQJLrVoyZIlXi8NAAAAyYDAABeVWmSFyN99952KFy+u3377zc1AePPNN3X69GmvlwcAAIBLQGCAi1a+fHmtWrVKbdq00cmTJ/Xss8+6Dkb//POP10sDAABAEhEYIEmyZ8+uL774QoMHD1aGDBk0ZcoUlS1bVgsXLvR6aQAAAEgCAgNcUmrRQw89pGXLlqlkyZL6448/VKdOHfXt25fUIgAAgCBDYIBLduONN2rlypW6++67derUKfXo0UNNmjTRrl27vF4aAAAAEonAAMkia9as+uyzzzRs2DBlypRJM2fOdKlF8+bN83ppAAAASAQCAyRratH999/vpiOXKlVKO3fuVL169dS7d2+3kwAAAIDARWCAZFe6dGlXd3Dfffe5WoNevXqpYcOG+vPPP71eGgAAABJAYIAUkSVLFn388ccuvci+nzNnjsqUKaNZs2Z5vTQAAADEg8AAKapDhw5asWKFbrjhBleMbDsHL730kpt/AAAAgMBBYIAUd+2112rp0qWutWl0dLT+7//+z9UeWHtTAAAABAYCA6QK61Rkw9BsKJp1MFqwYIHrWjR9+nSvlwYAAAACA6S2tm3batWqVSpXrpz++ecfN+/g+eef14kTJ7xeGgAAQFgjMECqK1GihBYvXqwuXbq4+2+88YabmPzbb795vTQAAICwRWAAT2TMmFHvv/++vvzyS2XPnt0FCraLMGnSJK+XBgAAEJYIDOCpO++8U6tXr1bFihW1Z88etWjRQt27d1dUVJTXSwMAAAgrBAbw3FVXXaVFixapW7du7v7bb7+tm266Sb/++qvXSwMAAAgbBAYICOnTp9c777yj8ePHK2fOnG5ysqUWjRs3zuulAQAAhAUCAwSUW2+9VWvWrFHVqlW1f/9+3X777eratauOHz/u9dIAAABCGoEBAk6RIkXcnINnnnnG3R84cKBq1KihzZs3e700AACAkEVggIAUGRmpfv36afLkycqTJ49WrlzpUovGjBnj9dIAAABCEoEBAlrTpk1dalHNmjV18OBBtWnTRp07d9bRo0e9XhoAAEBIITBAwCtUqJDmzp2rF154wd3/8MMPXQ3CL7/84vXSAAAAQgaBAYJCunTp9Nprr2n69Om67LLL9P3336tChQoaOXKk10sDAAAICQQGCCqNGjVyqUV16tTR4cOHdffdd6tTp046cuSI10sDAAAIagQGCDqXX365Zs2apZdfflkREREaNmyYKleurJ9++snrpQEAAAQtAgMEpbRp0+o///mPCxAKFCigdevWqVKlSho+fLjXSwMAAAhKBAYIajfffLNLLapfv75LJ7rvvvvUsWNHHTp0yOulAQAABBUCAwS9/Pnzu6LkPn36KE2aNPrss8/c7sEPP/zg9dIAAACCBoEBQia16MUXX3RtTa0GYf369a7uYOjQoYqOjvZ6eQAAAAGPwAAhpVatWi61qHHjxjp27JgefPBB3XXXXW44GgAAABJGYICQY3MOpkyZotdff93tJHzxxRdu5oEFDAAAAIgfgQFCktUaPPfcc1qwYIEKFy6sjRs3umnJgwYNIrUIAAAgHgQGCGnVq1fX6tWr1bx5cx0/flxdunRR69attX//fq+XBgAAEFAIDBDy8uTJowkTJuitt95SunTp9NVXX6l8+fJasWKF10sDAAAIGAQGCAs2Ifmpp57SwoULVaRIEW3ZssXtJgwYMIDUIgAAAAIDhJsqVaq41KLbbrtNJ06cULdu3XT77bdr7969Xi8NAADAUwQGCDu5cuXS119/rffee0/p06fX+PHjVa5cOS1dutTrpQEAAHiGwABhm1r0+OOPa/Hixbrqqqu0bds21axZ09UhnD592uvlAQAApDoCA4Q1m2+watUq16no5MmTevrpp9WiRQvt3r3b66UBAACkKgIDhL0cOXJo1KhR+uCDD5QhQwY3HK1s2bJatGiR10sDAABINQQGgD+16JFHHtGSJUtUokQJ/f7776pdu7abnkxqEQAACAcEBsAZbKdg5cqVat++vU6dOqUXXnhBt9xyi3bt2uX10gAAAFIUgQFwlmzZsmnEiBEaOnSoMmXKpBkzZriAYf78+V4vDQAAIMUQGAAJpBY98MADWrZsma677jrt3LlTN998s3r37u12EgAAAEINgQFwHtdff72WL1+ue++919Ua9OrVS40aNdKff/7p9dIAAACSFYEBcAFZsmTRJ598ok8//VSZM2fW7NmzXWqRfQUAAAgVBAZAIt1zzz1asWKF20X466+/1KBBA7388stu/gEAAECwIzAALoLVG1jdwYMPPqjo6Gi9+uqrqlevnnbs2OH10gAAAC4JgQFwkaxT0UcffaSRI0cqa9asWrBggcqUKaPp06d7vTQAAIAkIzAAkshmHdjMAwsK/vnnHzVp0sTNPSC1CAAABCMCA+ASlCxZ0k1LfvTRR919m5Rcp04dbd++3eulAQAAXBQCA+ASZcyYUf/97381ZswYZc+eXYsWLXJdiyZPnuz10gAAABKNwABIJq1atdKqVatUoUIF7dmzR82bN9fTTz+tqKgor5cGAABwQQQGQDIqXry42zHo2rWru//WW2+pVq1a2rp1q9dLAwAAOC8CAyCZZciQQQMGDNC4ceOUM2dOLV26VOXKldP48eO9XhoAAECCCAyAFNKyZUutXr1aVapU0b59+3TbbbfpiSee0PHjx71eGgAAwDkIDIAUVLRoUTfnoHv37u7+e++9pxo1amjz5s1eLw0AACAOAgMghaVPn179+/fXpEmTlDt3bjf7oHz58vryyy+9XhoAAEAsAgMglTRr1kxr1qxxOwYHDhxQ69at3fyDY8eOeb00AAAAAgMgNRUuXFhz5851E5LNBx98oKpVq2rDhg1eLw0AAIQ5AgMglUVGRuq1117T9OnTlTdvXq1du9bNPvj888+9XhoAAAhjBAaARxo1auSCgtq1a+vQoUO666679OCDD+rIkSNeLw0AAIQhAgPAQ5dffrlmzZqlnj17KiIiQkOHDnXtTX/++WevlwYAAMIMgQHgsXTp0ql3796aOXOm8ufPrx9//FEVK1bUp59+6vXSAABAGCEwAAJE/fr1XdeievXquXSie++9190OHz7s9dIAAEAYIDAAAkiBAgU0Y8YMt4OQJk0at2tguwe2iwAAAJCSCAyAAJM2bVpXczBnzhwVLFhQ69evV6VKlVz9QXR0tNfLAwAAIYrAAAhQ1q3IUouse5ENQbOORXfffbcOHjzo9dIAAEAIIjAAAli+fPk0depU9e3b1+0k2KwDSy2ygAEAACA5ERgAAc5qDZ5//nnNmzdPhQoVclOSbVqyTU0mtQgAACQXAgMgSNSsWdPtFDRt2lTHjx/Xo48+qrZt22r//v1eLw0AkIrsktCfkjZIWitpqaRlkn6QtFHSHq8XiKBFYAAEkTx58mjixInq37+/m38wZswYlS9fXitXrvR6aQCAFAoCtksaJ+klSU0k5ZdUUNI1kspKqiqpiqQbJZW0fyskFZV0h6TXJM2UtNvrXwRBgcAACMLUou7du+vbb79VkSJFtGXLFlWvXl0DBw4ktQgAQsQqSQ/6g4ArJd0u6f8kTZf0dyKev03SWEkvSmokKa+kYpJ6+I8B8SEwAIKU1RmsXr1aLVu2VFRUlLp27ao77rhDe/fu9XppAIAkOCbpf5KqSaogaWgig4DE2iqprz9AaOEPMk4n4+sj+BEYAEEsV65cGjt2rAYMGKDIyEiNGzfOpRYtXWoZpwCAYLBT0nOSCkm6R9KSFH4/21ue5E9LKiGpv6QjKfyeCA4EBkCQi4iIcLsFixcv1lVXXaWtW7e6QuW3336b1CIACGDR/h2CUpL6eVQHsEXSM/76hG89eH8EFgIDIETYfINVq1bpzjvv1MmTJ10dQosWLbR7NyVnABBodvjTeWyHYJ/Xi5G02QZrSurG7kFYIzAAQkiOHDlcp6JBgwYpQ4YMmjx5ssqVK+d2EwAAgbNLUFrSZAXe2gawexDWCAyAEEwt6ty5s5YsWaISJUpo+/btqlWrlt544w2dPk2ZGQB45Yi/hWig7BJcaPfAOhiRkBpeCAyAEFW2bFk336Bdu3Y6deqUm55sw9H+/js5e1wAABLDAoGG/nkEwSDa38HIWqae8noxSDUEBkAIy5Ytm0aOHKkhQ4YoY8aMmj59ugsYFixY4PXSACBs2OWYupIWKfgMk9RO0gmvF4JUQWAAhEFqUadOnbRs2TJde+212rFjh+rWras+ffq4nQQAQMrZL6mxpDUKXl9K6sjOQVggMADCxA033KDly5frnnvucbUGPXv2VOPGjfXXX395vTQACNmagub+KcbB7gtJj1FzEPIIDIAwkjVrVn366af65JNPlDlzZs2aNUtlypTRnDlzvF4aAISch0Osu8+Hkt71ehFIUQQGQBi699573e5B6dKl3Y5B/fr11atXL1KLACCZjJc0QqHHOhX94vUikGIIDIAwVapUKVd38MADD7gJyb1793YBgtUgAACSzsZKPqLQdEzSfdQbhCwCAyCMWTrR0KFDNWLECGXJkkXz5s1zXYtmzpzp9dIAIGg9LimUq7e+8w9CQ+ghMACgu+66S6tWrXL1BjbnoFGjRurRo4dOnjzp9dIAIKiM8xfqhroXSSkKSQQGAJySJUvqu+++0yOP+DbA+/bt69qa/v77714vDQCCwt4QTiGKL6XofroUhRwCAwCxMmXKpA8++ECjR492w9EWLlzoUoumTJni9dIAICiGge1S+Fgsab7Xi0CyIjAAcI7WrVu71KLy5ctr9+7datasmZ555hmdOMHsSwCIz2lJHyj8DPJ6AUhWBAYA4nX11Vdr8eLFevxxK6OT+vfvr1q1amnbtm1eLw0AAo61bNii8KypoJdd6CAwAJCgDBky6L333tPXX3+tHDlyaMmSJS61aMKECV4vDQACSrheObcWFUO9XgSSDYEBgAu6/fbbtXr1alWuXFn79u1Ty5Yt1a1bN0VFRXm9NADw3FZJkxW+Bksi0TQ0EBgASJRixYrp22+/1VNPPeXuDxgwQDVq1NCWLeG4eQ4A/xoS5t15LJVokteLQLIgMACQaOnTp9dbb72liRMnKleuXFqxYoXKlSunr776yuulAYBnpnu9gAAww+sFIFkQGAC4aM2bN9eaNWtUvXp1HThwQK1atVKXLl107Jh1tgaA8HFc0g9eLyIArPR6AUgWBAYAkuTKK6/UvHnz9Nxzz7n7gwYNcoHCxo0bvV4aAKSaH8mvdyw4ouos+BEYAEiyyMhIvf7665o6dary5s3rCpRt9sGoUaO8XhoApAqulPtE+YMkBDcCAwCXrEmTJi616KabbtKhQ4fUrl07Pfzwwzp69KjXSwOAFLUiMT+0dasUESHVqSMdPixZE4fChW3cvFS+vDTpjNLdL7+UqlSRsmSR8ueXunaVzv67dM0a6dlnpQoVpMsus97S0lVXSY8+Ku3Ycf73P3DA9/7FitnVHalbN+mxx3zHP/oo4d/hmmukNGmk8zScSNRngYBGYAAgWVxxxRWaM2eOXnrpJUVEROijjz5SlSpVtH79eq+XBgCBsWNgLZ7r1ZNGjpSqVvXd1q6VbrtNmjVLeucdqX17KVs2qVEj6dQpaeBAqVOnuK/z+uu+nzU1a0q33CJFR0sffCBVrBh/cGAswKhdWxo+XCpbVmrRQsqVS3r4Yd/xIdZfKR7z50sbNvjWbgFIcnwWCEgR0dH2X9LFK126tPu6bt265F4TgCD3zTff6O6779auXbuUJUsWffDBB+rQoYPXywKAZJdJ0gXbLtgVe7tCb26+WZo40bcjYOwk/b77bNy8tHu3NHOm7+Te2Al+uXLSrl3S5s3/npTPnSuVKuXbUYhx+rTUp4/Uq5fv9T7+OP73r1ZNmjpVypkz7hpr1JAWL5ZWr/YFDWe6+25fMDN6tNS6dYK/Zg1JCy/0WSBFXer5OTsGAJJdgwYNXGpR3bp1dfjwYd1zzz26//773fcAECpOJSYoOJOl4thV/ZigwNxzj5Q3r7Rpk9Sly79Bgbn8cumuu3zfL1jw7+N168YNCmJe++WXbfvWF3gk5L33zg0KzCOPxL9rsHev9PXXvpSlli3P++vxN3zwIzAAkCIKFizodg7+85//KE2aNPrkk0/c5GR2GQGEiotu0Fy0qFSy5Lkn9EWK+L5v2PDc58TsEuzcGfdx21345BOpe3fpgQeke+/13U6c8B3bs+fc1ypYMG7gcaZWraQ8eXw7A0eO/Pv4iBGStaLu2NGG2Zz316NhdfAjMACQYtKmTauXX35Zs2fPVoECBfTTTz+pUqVK+vjjj5XELEYACBgXfSJsV/PjkzVrwsdjjh23iQl+X3zhCzLuv196+21f2tCnn/pulnZkDh4897WuvDLhtWXM6Dv537/fVwAdY+hQ39ez6xziQWAQ/AgMAKS4OnXqaO3atWrYsKHrVPTAAw+49CLrYAQAwSrdxT7Bdgcu5bjZts23M2CFzO++K9nsGLvCbxdb7GY1BCa+iy928n8+VoRs3Yli0omWLZO+/16qVcvXlSi5Pw8EHAIDAKkiX758mjZtmv7v//7PpRaNGDFCFSpU0Pf2jw4ABKELnGanDCsctqDA2pg+8YSvaNnansY4TzvRC7I0J6tfWLRI+vnnfwOEhx4K3M8DyYrAAECqsYCgR48ebmKytTfdsGGDqzsYPHgwqUUAgo5l3Eek9ptaMbApVOjcY1ag/Ndfl/b6MUXIlqJkwyqtnekddyTqqWeEJwhSBAYAUp0NQrOuRbfccouOHz+uRx55xA1FO2CDdwAgSFhQUCC13zSmeNmKgs/s9PbHH/+e1F8K6zxUoICvtsDSPa3V9IVSkPzO6pOEIERgAMATefPm1aRJk/Tmm28qXbp0Gj16tEstWrVqlddLA4BEK5/ab2hDyaxX/YoVvjSiO++UmjXzBQx2db969Ut7fZuGbEXNMRKZRmQqXNo7IwAQGADwNLXo6aef1oIFC3TllVdq06ZNqlatmt5//31SiwAEhQSaf6Ycaxn67bdS586+K/mTJ/vqAR5/3KZL+k7sL5UNYTNWyOwfmBWQnwWSHZOPAQSEPXv26L777tNE/2CeO+64Q0OHDlXO+AbxAECAmGQX8RVirDvRRx/55iRYB6RE2mGjElJ0YbgQJh8DCAm5c+fW+PHj9e677yoyMlJff/21ypcvr+XLl3u9NAAIn/QZa4f6v//5pjG3aZPop1lAQFAQ/AgMAASMiIgIPfHEE1q0aJGKFSumX3/9VTVq1HDBAqlFAALR5ZIuO3lSQe/NN32FxlWqSEePSi+9FLcNargFSGGKwABAwLHpyFaEbOlEJ06c0JNPPqmWLVu6dCMA8JpdqPjxxx/Vp08f1zThb8vzD3ZTpvg6HaVNK/Xq5ZuTcBEqp9jCkJqoMQAQsOyvp0GDBumpp55SVFSUChcu7LoXWYEyAKSm06dPa+nSpRo3bpy7WbOEGBF33KHor75SOLdt3SypmNcLgagxABDSqUVdunTRkiVLdPXVV2v79u1uBkK/fv3cP9IAkJLsgsSMGTPcrBUbyli9enXXYtmCggwZMqhZs2YaNmyY/vjgA5dSFK5uISgIGewYAAgKNvzs4Ycf1iibxCmpSZMm+uyzz9w8BABILocOHdK0adNcM4QpU6Zo//79sceyZ8/uggFLbWzcuLGyZcsWe6y3pF4KT1P8wQG8d6nn5wQGAIKG/XU1ZMgQV6B87NgxdwXviy++cLsIAJBUf//9txu4aClC33zzjZvIHqNAgQK69dZbddttt6lu3bpKb3MEEmjVWURSCJQhXxTbKdgoKa3XC0GynJ+nS9KzAMCj1KKHHnpIVatWVevWrfXLL7+oTp066t27t1544QU3MA0AEmPr1q1uV8CCgYULF8ZJT7TURQsE7FalSpVE/d1iqUS3SfpS4eURgoKQwo4BgKDd7n/00Uf1P+u3LalBgwbu+/z583u9NAAB3EnIAgELCFavXh3nuM1NiQkGSpUq5S5EXKwFkmorfGSUtF0SCZ2Bgx0DAGEpa9as+vTTT93WvhUo2/Z/2bJl9fnnn7vHAMB2Ab777rvYnYHNm613jo/tAlgaogUCVjNQpIglAl2aWpJsJNhohYdXCQpCDjsGAIKe/T1kqUU//fST+8f+5Zdf1ksvvaS01o8bQFix+oC5c+e6QGDChAn666+/Yo9ZJ6GGDRu6YKB58+Yp0rzgHztHkrRLoc2aRn9LGlHAofgYACQdOXJEjz/+uD7++GN333YNRo4cqYIFC3q9NAAp7ODBg66TkAUDU6dOdV3MYuTIkcN1ErJgoFGjRm63MaWNlXSHQlcGSWslXeP1QnAOUokAQFLmzJldP3ErRu7cubO7YmipRSNGjHD1BwBCy65duzRx4kQXDMyaNcvNHIhhFwQsPchu9ndCQp2EUsrtktpK8jVXDj19CApCFjsGAELO+vXr1aZNG33//feugLBHjx565ZVXlC4d10KAYPbrr7/G1gssWrQoTiehEiVKxBYPV65c2fMuZaGaUkQKUWAjlQgA4nH06FE9+eSTGjx4sLtvRYZWmFyoUCGvlwYgkewU5YcffnCBgN3WrrUEln9VqFAhNhi47rrrktRJKCXNkNQshGYbWEXGd9bO1euFIEEEBgBwHjYp2WYfWA5ynjx5XEtTm5oMIDCdOnXKdRKKaSu6ZcuW2GPWUKBWrVouELChY1deeaUCnaUTtbcgR8Etu6S51tbV64XgvAgMAOACNm3a5LoWxfQtf/bZZ9WnTx9FRkZ6vTQA/k5Cs2fPdsGA1Q1Y/UCMjBkzuqJhqxewTkIW4AebjyQ9rOCVyb/7wYz5wEdgAACJcOzYMT3zzDN6//333f1q1aq53YRguOIIhCLrHGQdhGxXwL7arl6MnDlzxukklCVLFgU765f2oM1WUHDJIWmKpBpeLwSJQmAAABfh66+/1gMPPKD9+/crV65cGj58uFq0aOH1soCwYDMFYjoJ2Q7BmZ2ELr/8crcrYMFA7dq1Q3JH72tJ7SSdUHDI598pKOv1QpBoBAYAcJEsZ7lt27Zavny5u29Fyq+//nqqtzQEwuX/bzHFw4sXL3YFxTFKliwZWzxcqVIlzzsJpYalku617mkKbPUlDZV06fOgkZoIDAAgCexK5XPPPad3333X3bf2hpZaVKxYMa+XBgQ1O62w7kExwYB1FTpTxYoV43QSCkfHJPWS1D8AU4uy+tf1kJ0ker0YXDQCAwC4BJbWcO+992rv3r1uQqpNTr79dhtPBOBiOgnZXIGYGQNbt26N00nIUoNiOgkVLlzY07UGkiWS7gug3QN2CYIfgQEAXKJt27a51KIlS+yfaemxxx5T//79lSFDBq+XBgR0Qf+ZnYT+/vvv2GOZMmWK7SRkRcTB2EkoNXcPXpFke5fHPVpDLkl92SUICQQGAJAMTpw4oZdeekn9+vVz98uXL6/Ro0fr6qsZ5QPEsKJ96yBkwcC0adN06NCh2GNWzB/TSahhw4Yh0UkotSclfyLpA5vwnErvWUnSo5La+FuSIvgRGABAMrKTnnvuuUe7d+9WtmzZNHToUDcDAQhXf/75pyZMmODShGyHwILoGFdccUVsJyEbPBaKnYRS22l/J6D/2t9HKTAYLaO/M1Jnf2CA0FKawAAAktfvv/+udu3aaeHChe7+ww8/rHfeecelRwDhMhQwpl7AphCfeapw7bXXukDAAgIrJA6HTkJesZ2DyZJW+m8/JaFY2RIiy1jRt/92q6TcKbReeI/AAABSwMmTJ9WrVy/17dvXnRTdeOONGjNmjK655hqvlwYkO/tvfM2aNbGdhH788cc4x61rV8zOgAUG8MYRSWv9QYJ9tZFwR/03C88y+XcErKKjvKQKkkpJYh8nfJQmMACAlDNz5kzdfffdrrDScqY//PBDdx8IhU5CtitmgYDtDlgR/pmdhOrUqRPbSahQoUKerhVA6pyfp0vSswAgTFgRpfVkv+uuuzR37lx16NDBfR04cKAyZ87s9fKAi+4k9M0337hgYNKkSfrnHyt59bFUucaNG7tgoGnTpsqdm4QTINywYwAAiby6+uqrr6p3794u7cL+DrTUolKlbKMeCFz79u2L00no8OHDscfs5L958+YuTciCYIJdILiRSgQAqWjOnDlu98A6tdgV1kGDBrkBaUAg2blzp+skZMGA7XCd2UnI0oLO7CSULh3JA0CoKE0qEQCknptvvtkVaVpKkaVk3Hfffe7E67///a+yZs3q9fIQxjZu3BhbPBwzrC/Gdddd5wIBu1WoUEEREYyxAnAudgwAIAlOnz7tOha9/PLL7nvr1GKpRTfccIPXS0OYsH++V61aFVs8fPa/x1WqVIltK0o3LSA8lCaVCAC8s2DBAjfzYMeOHcqYMaPee+89derUiSuySLE2umd2Evrtt99ij1lKUN26dV0gYJ2EbPgYgPBSmlQiAPCO5WhbalHHjh1dYedDDz3k6hAGDx6s7NmzK1DttaJU61Lj74Ee4e9/bn3Qc0nK4fUCEevo0aNxOgnZVO4YVix8ZiehXLnsTw8AkoYdAwBIBpZO1L9/f/Xo0cN1MLr66qtdalG5cuW8Xpp2nTE51W4rbLrzBZ5TzD8ltYL/ZsOSaF6Zup2EJk+e7IKB6dOn68gRG231byehFi1auGCgQYMGTOQGEItUIgAIIIsXL1bbtm21fft2pU+fXu+88446d+6c6qlF6yV9IGmcpO3J9JpXSWot6WFJRZPpNfEvS0ez9CC7WUG7pQ3FKFy4cGzxcM2aNekkBCBeBAYAEGD27NnjWpha2oe58847NXToUOXIkbIJOtaQcqKkQdZWNQXfx0KcppIeldRIUpoUfK9Qt2HDhthOQkuXLj3n39mYtqLly5enbgXABREYAEAAsr9a3333XT333HOuh3yxYsVcalHFipagk7z2SBoo6SO76qzUZbsIj/iDhCyp/N7B+t/FypUrY4OBn3/+Oc7xatWqxQYDJUqU8GydAIITgQEABLBly5apTZs22rp1qyIjI/Xmm2+qa9euyXb1d6I/tedPecsChI8l1fZ4HYHIUoKse1VMJ6Hff/+3wsNSgmw2hgUC1kmoYMGCnq4VQHAjMACAICgkfeCBBzR27Fh3364If/zxx5fUQcZ2CZ6QNEKB5TFJr7N74IqFZ86c6QIBSymz9LIYWbJkUZMmTVwwcMsttyhnzpyerhVA6ChNYAAAgc/+qrXpyN27d1dUVJSKFCmiUaNGqWrVqkG7S5CQcN092Lt3b2wnoRkzZsTpJJQ3b17XSciCwvr169NJCECKIDAAgCBi+eWWWrR582aXRmLTk5966imlSXPhEt4ofz7/JwoOT0t6I8SLk//444/YTkLz5s2L00nIgr+YeoEaNWrQSQhAiiMwAIAgc+DAAT344IOuGNnYYKpPP/1UefLkSfA5du35DknTFVzu9gcyoXRKvH79+th6AashOdP1118f21a0bNmydBICkKqYfAwAQcYmIlsakRWdPvHEE5oyZYo7ifziiy9cj/qzHfC3B12o4GM1EAcljZaUQcHJrp8tX77cBQIWEFhgEMNO/K2TkAUCtjtgg+0AIFixYwAAHlq7dq1at27t+tmnTZtWr776qmtxGpNadFRSY0kLFNxst2NUEF2NshazZ3YSspShGNZdql69ei4YsLqBAgUKeLpWAIhBKhEABLlDhw656cgjRvh6DDVq1EifffaZcuXLp9skTVFouFfSsACuObBiYSsatmDAioitmDhG1qxZXQch2xWwryk9rA4AkoJUIgAIcnbSaYFA3bp19dhjj7mTU0stqrZsmaYUKqRQMVySjezqocBhbUTP7CR09Kjt0fhcdtllbkfAdgZshyBjxoyerhUAUhqBAQAEAMtVv//++1WlShW1atVKP+fPr7EhFBTEeEVSc0k3eLgGGzAWUy8wf/58nTp1KvZY0aJFY4uHq1ev7tK7ACBcEBgAQIBtA89dvlxXHzmiQwo9J/wpRUssVz8V3/fnn392gYDdVqxYEefYjTfeGFs8XKZMGToJAQhbBAYAEGBezZJFh7KE7uzgVZL6SXoxkT//zz//6IcffnCpVol1+vTpOJ2Efvnll9hjduJvcwUsELBb8eLFk/BbAEDoITAAgAAyV9J/Ffr+I6lFIlKKrFuT5fdb+s/UqVPVpEmT83YSstSgmE5CO3bsiD2WPn36OJ2E8ufPn4y/DQCEBgIDAAgQNsTsAYWHmJQiGw+WUBb/Tz/95E7m//zzT3ffZj+cHRgcPnw4Tiehffv2xR7Lli2b6yBkwYA9z+ZHAAASRmAAAAHic0m/Knys8k9ytuFt8c13qF+/vksjijFp0iS3K2CTo+17CwZmzpypY8eOxf5Mvnz5dOutt7pgwAbIZcgQrGPVACD1ERgAQACIDpMUorMNiicwsOLghg0bxpkjYOx+pUqV9OOPP8bpJHTVVVfFFg/bFGI6CQFA0hAYAEAAWCppjcLPNElb7OTef/+7775T48aN3a6AEthJMNY9KKat6A033EAnIQBIBgQGABAgV87DdadksKQ3rPB67lw1bdo0zpCxs9nE4VWrVrldAgBA8grUyfQAEDYsi360wtcwSV9NnuwKjc8XFJj9+/dr9+7dqbY2AAgnBAYA4LFPJEUpfNlp/sjjxxUdbfsHF2ZFxwCA5EcqEQB4bOqlvsC8eZIN/+rYURo+XMEo8x13uA5Dv/32mzZt2qT169e7acX2/ZmFxmbhwoWerRMAQhmBAQB46LS/bWe4W2m7Bg0anPN4VFRUnEBh8+bNuv322z1ZIwCEOgIDAPDQZknx998JLxskHbShZGc9bhOLS5Uq5W4AgJRFjQEAeHylHL7uRKu9XgQAhDkCAwAIxMBgyhTp/vul666TsmeXsmSx5v3Sa69Jx48n7sX37ZMGDpQaNZKKFJFsCnCePFLjxtI338T/nDp1JJsJsHWr9PnnUtWqUrZsUs6cvuOvvOI7brUMK1dKTZr4juXOLbVuLf3+u+/nDh+Wnn1WKlpUyphRuv566auvkvZZAABSBYEBAHhoRUIHHnhA+vpr3wm3nXzfdJO0fbv04ovSLbdIZxXkxmvJEqlrV2nDBumaa6TbbvN9nTnTFyx8/HHCz+3bV+rQwXJ5pGbNfCf2Z1q6VKpRQ/r7b99rWcDx5ZdSvXrWU9RXDP3pp1KlSlK1atJPP/kChxkzLv6zAACkCmoMAMBD6xI6MHiw1LChlCnTv48dPCi1by9NniyNHCndc8/5X9yCgO++8131P9Pq1dLNN0tPPuk7Wc+a9dznfvaZNGeOVLt2/K/94YfSBx9Ijzziu3/ihC9gmTVLql5dKlBA2rLFt9Nhhg2TOnXy7XhYIHExnwUAIFWwYwAAHjqc0IFbb40bFBhL6XnnHd/3EyZc+MWLFTs3KDDlykldukgHDti44YR3LBIKCkzNmv8GBSYyUnr8cd/369f7goaYoMDce6+UN68vULEg4mI+CwBAqmDHAAA8LLg975zfjRulqVOlTZt8OfunT0sxQ8DsWGJYytHs2dLixdLOnf/WJ8Q8P6HXadHi/K9ruxlnu+oq31erKyhZMu6xtGl9dQ5Wl/DPP1LBguc8/fwzjwEAKY3AAAA8YtfN4531ayf/Tz/t2x1IaBqwpRVdiBUCW33A2rUJ/0xCr3Plled/7SuuOPexmJSk+I6deTyB4mkCAwDwFqlEAOCRtAkdGD1aevttqVAhXyefP/6wSV++ICHmpDqhgOFMltNvQcEdd/iKha1Lke0g2HOthuF8r2OdhM4nTZqkHTsPrlQBgLf4exgAPAwMIv07B3GMG+f7ann6TZvGPWYFvYlhqUfWkjR/fl+gYak8SXmdVHRWRQUAIJWxYwAAHor3uvzevb6vtmNwtjFjEvfC1jLUahIsl//soMCKf2OCjwBygT0KAEAKIzAAAA8ViO/BmMLdjz6Km+rz7bfSm28m7oXz5ZNy5JB+/FFatOjfxy2V6LnnfLMNAkx+rxcAAGGOwAAAPFQ+vgdtKJm1+hw0yDdYrF07qVYtX/vQM1uEnk+6dL7JwydP+p5nXYTatpWuvto3g8DalQaYCl4vAADCHIEBAATaybDtGKxYITVv7mvtOXGidOiQr2A4sTsGpkcP3/ThG2/07RrY8LEyZXwTkStWVKAhMAAAb0VERyemtcW5Spcu7b6uW8esSgBIKhsvdrPXiwgQ621Ys9eLAIAgdqnn5+wYAEAKsGsuR48eTVoqURjKJqmE14sAgDBHu1IASAEtWrTQ5MmTVbhwYV177bW67rrrYr/aLV++fIqIiFAOSVdL2qTwVo4rVQDgOQIDAEgBu3btcl+3b9/ubt/YTIEz5MiRQ7lz59Zll12mMvPna9OFBoqFuMpeLwAAwAUaAEipHYPz2b9/v3799VctW7ZM2caPV7jr4PUCAAAEBgCQEm677bZE/dyNN96oD26/XcUVvmrY5+D1IgAABAYAkJxOnjypOXPmaNCgQUpnswTOo02bNlqxYoUypk+vzgpfj3q9AACAQ40BAFwi6z40c+ZMjRs3TpMmTdKePXsu+Jx77rlHH3/8sdKmTevu3yvpJUnHFF4uk3SH14sAADgEBgCQBPv27XNdhywYmD59uo4cORJ7LE+ePK7GwLoPPWvTh8/SqVMnDR48WGnS/Ltpm0dSW0nDFV46Scrg9SIAAA6BAQAk0s6dOzV+/HgXDMydO9elDcWwtqRWV2C3mjVrujSi06dP65133nHPi/HYY49pwIABcYKCGF3CLDCwvZKHvF4EACAWgQEAnMfGjRtdIGC3JUuWxDlWqlSp2GCgfPnybi7Bmezkv2XLlvrggw/c/e7du+vNN9885+diVJR0n6RPFB6ek1TU60UAAGIRGADAWROLV69eHRsMnD1WvkqVKrHBQMmSJS/4el26dNGsWbNcTcGLL76YYFAQ421JMyX9odB2vaSXvV4EACAOAgMAYe/UqVNauHChCwQsVWjbtm2xxywlqE6dOi4QuPXWW3XFFVdc1GuXLl1aGzZsSPTP55Q0RNItCl1po6P1SUQEtQUAEGAIDACEpWPHjrkr+RYMTJw4Uf/880/ssUyZMqlx48YuGGjWrJly5cqVqmtrEuIpRRkHDND+G26Q6tXzeikAgDMQGAAIGzZteOrUqS4YmDZtmg4dOhR7zE7+mzdv7oKBhg0bKnPmzJ6uNVRTijJs3KjDzz2nBidO6KWXXtLLL798wXkPAIDUwd/GAELaX3/9pQkTJrhgYPbs2Tpx4kTsMUsLsuJgCwZq1aqlyMhIBQpLKRohqZGkKIWG7JJmFCqkYffco6FDh+rVV1/V/Pnz9cUXX+jyyy/3enkAEPYIDACEnC1btsQWDy9evNgVFMe45pprYouHK1asGG/b0EBRR9IoSXdKOq3gllHSZElVM2VS1SFDVLduXT388MNasGCBypQpo//9738ufQsA4J2I6DP/xbzIgjpzdscOAEht9tfY999/HxsM2PdnsgAgJhiwoWPB5lN/zUGS/rIOALYPMz6egmorym7durXWrl3r7j///PNuF4HUIgBImks9PycwABC0nYS+++672E5CtksQI23atC41yAIBSxWy4WPBztKK7rXfW8Elk6Rx/pSohIrAbb7DoEGD3P0aNWq41KJQ+DMDgNRGYAAgbBw/flxz5sxxwYDVDezatSv2WMaMGV3RsAUDVkScJ08ehZqJklrb56DgqSmYIqlmIn72yy+/VKdOnXTgwAHlzp1bn376qesIBQBIPAIDACHt4MGDroOQBQPWUchOHGPkyJHDnTxaMGD56VmyZFGoWyGpo6SfFNiqSBou6dqLeM7mzZvVpk0brVy50t1/6qmn1LdvX6VPnz7F1gkAoaQ0gQGAUPP333+72QIWDNisAdspiFGwYEE3aMyCARs8Fo4njfZp/EfSGwFYlGxDy161k3pL6UrC8+3P+tlnn9V7770XO2l61KhRKlq0aLKvFQBCTWkCAwChwKYNxxQP2xTi06f/PeW9+uqrY4uH7UQxkDsJpabl/rqDnwJol8CGsiVHebfVjdx3333at2+fcubMqU8++cTViwAAEkZgACAo2V899vdHTDCwevXqOMfLlSsXGwzY3zcRERGerTUYdg/e9rD2IKuknpK6J3GXICFbt25V27ZttXTpUne/a9eu6tevnzJksH0JAMDZCAwABA3bBbCTvJhgYNOmTbHHbBegZs2asZ2ESB25OP/4r9Z/IOnXVHrPUpK6SLrbX2icEqKiotSjRw+99dZb7n6FChU0evRoFS9ePIXeEQCCF4EBgIBmJ3bz5s2L7SS0c+fO2GN25bdBgwYuEGjRooUuu+wyT9caCiwBa4akQf6OQMk9+8AmDNwu6VFJtewfEaWOyZMnq2PHjtqzZ4+yZ8/uJie3atUqld4dAIIDgQGAgHP48GFNnz7dBQN2Qrd///7YY3ZS17Rp09hOQtmyZfN0raHsV3+L05X+2/okFCtbatD1dqXeBsVJsiz/gvLG9u3b1a5dOy1atMjd79y5s95++23XqhYAIAIDAIFh9+7dmjRpkgsGZs6c6QZXxcifP39sJ6G6deuSI+6Rw5LW+Fue2qxhC9fsT+mo/8p/Rv9AslxW4+EPBm70PxYoTpw4oV69erk2pqZMmTIaM2aMSpYs6fXSAMBzBAYAPGNXcK17jAUDCxYscNOIYxQrVky33367CwaqVq3qphEDyWXGjBm6++679c8//yhr1qwaPHiw2rdv7/WyAMBTl3p+bumiAJBoP//8c2zx8IoVdu35XzfeeGNsJyH7nk5CSCmNGjXS2rVrXTAwf/583XXXXZo7d64GDBigzJkze708AAhK7BgAuGAnIQsAYoKBX375JfaYnfjXqFEjtpPQVVdd5elaEX5OnjypV1991d3sn7Prr7/epRZdd11yTFMAgOBCKhGAFMnjttQgCwQsVeiPP/6IPRYZGan69eu7YMA6CVn9AOC12bNnu12Dv/76y+0YDBo0yHUxAoBwUppUIgDJ4ciRI65o2IIBKyLeu3dv7DHL4b7llltcMGBfrbMQEEjq1aunNWvWuLoDCxLuvfdel1r03//+V1myZPF6eQAQFNgxAMKYnfxbO1ELBqy96NGj1p/GJ2/evLGdhOyki5aQCAZWAP/aa6/plVdecWlw1157rb788kuXYgQAoa40qUQALsaOHTtiOwnZ4DHL0Y5RpEiR2OJhqx2gkxCClRUkW2Gy/fduQe3AgQP1wAMPUBAPIKSVJpUIXonpf25foySlP6MPOteWA8uGDRtii4eXLl0a55hdSY0JBsqWLcuJE0JC7dq1XWpRhw4dXGvTBx980KUWffjhhwzVA4AEsGOARNlxxvTUmNvO8/z85f7hSGfevJqWGo7s/9arVq2KDQZ++umnOMerVasW20moRIkSnq0TSGmWTtSvXz+99NJLLs3I/nu3rkUWBANAqClNKhFSwmlJcyQNkfTtBYKAxLJg4SZJD0mq65+0iuRjKUELFy6M7ST022+/xR5Lly6dbr75ZhcMWN1AwYKEaQgvixYtUtu2bfX777+7ydvvvPOOHnnkEXbIAISU0gQGSE7Wh+ZTSR9Y+kkKvs+1kjpLukdSzhR8n1B37NgxffPNNy4YmDhxonbv3h17zFo2NmnSxAUDTZs2Vc6cfNIIb/b/D+tWZAX3plWrVhoyZIhy5Mjh9dIAIFkQGCBZfC/pPUmf++sGUovNJ71LUlfLdU/F9w1m+/fv15QpU1wwMG3aNB0+fDj2WO7cud1sAQsGGjRooEyZrOIDQAz7J892C5577jm3y2ZD+UaPHq2KFSt6vTQAuGQEBrgkVjj8sqS3/OlDXkkj6VlJr0jK4OE6AtWff/6pCRMmuGBgzpw5bgBZjEKFCsUWD990000ubQjA+VkRfps2bbRt2zY3tK9///56/PHHSS0CENQIDJBkSyTdJ2m9AkcpScMlVfJ6IQFg8+bNscXD3333nbvSGeO6666LDQYqVKjAyQyQxDke999/v6vJMfb/p2HDhilXrlxeLw0AkoTAAEG7S5CQcN09sP8rrl27NjYY+OGHH+Icr1y5cmwwcM0113i2TiDU/n9nMw6efvpptxNnszwstahKlSpeLw0ALhqBAS7Kz5JuD7BdgvPtHoyVFMqnwNY+cfHixbHBwNatW2OP2XCxOnXqxHYSspQhACljxYoVLrVoy5YtLh3v9ddf11NPPcVuHICgQmCARFshqbF15lDwyCtphqTyCh3Hjx/X7NmzXSBgdQN///137DErFm7UqJELBpo1a+aKiQGkXmG/DUL78ssv3X37/+Dw4cOVJ08er5cGAIlCYIBEmS+puaSDCj7ZJU32z0AIVgcOHHAdhCwYmDp1qg4e/PdPwtqINm/e3AUDFhRYm1EA3rB/EgcPHqxu3bq5IN526kaNGqUaNWp4vTQAuCACAySqyLi+pH+bWgafrJJmW569gseuXbvcbAELBmbNmqWoqKjYY5dffrmbOmzBQO3atV1XFACBY82aNWrdurU2btzo0vr69OmjZ599VmnSWBUUAAQmAgNccD5BbUn7FPwsqWaB/benwGU1AjH1AjZp9fTpf8u7S5YsGVs8XKlSJU4wgABnO3s2Hfnzz23Ci9S4cWN99tlnuuyyy7xeGgDEi8AACfpH0o2Sdip0XOEPdgIl897+7/Pjjz/GBgN2lfFM1ko0JhiwFqMUMgLBxf4//vHHH+uxxx5zk8Ztt88CBdvpA4BAQ2CABLWVNFqh525J//Pw/W0XYMmSJbHBgM0biGG7ALVq1XKBgKUKXXnllR6uFEBysQsArVq10vr1693/z1955RX16NHDpRkBQKAgMEC8vpZ0p0LXBEktUvH9rD5g7ty5sZ2EbBJxjAwZMqhhw4YuGLAi4rx5rZcSgFBz+PBhdenSRZ9++qm7X69ePY0YMUIFChTwemkA4BAYIN4UIpsB8G8TzNBj/wyvS+GUokOHDmn69OkuGJgyZYprZRgje/bsrpWhBQOWd5w1q5VHAwgHFhg8+uijOnLkiPLnz6+RI0e6IAEAvEZggLBJIUqNlKJ//vkntpPQN99849oVxrCrgjZozIKBunXrKn369Mn87gCCxc8//+y6FlmKkdUO9ezZUy+//DKpRQA8RWCAOMb5JxuHi4n++QyX4rffftP48eNdMLBgwYI4nYSKFy8eWzxctWpVOgkBiGU7Bk888YSGDh3q7ltBshUmW4EyAHiBwACx7A+yjKQfFD7KSVpp/yFfxHPsP3m72hdTPLxypb3Cv8qWLRsbDFx//fV0EgJwXhYMPPzwwy790GqMrO7AhhUCQGojMECshUE+HfhSBrhVucDP2C7A8uXLY4OBDRs2xB6zE/+aNWvGdhIqVqxYiq8ZQGixv1MstWjt2rXu/vPPP69XX31V6dKl83ppAMJIaQIDxGgv6QuFn3usGDCex0+cOKH58+e7QMBShXbs2BF7zOoD6tev74KBFi1aKF++fKm6ZgChx+YcPPXUU/rggw/c/Ro1auiLL75Q4cKFvV4agDBRmsAA5i9J9k/PCYWfDJJ+l5TXn/M7Y8YMFwxMnjxZe/fujf25bNmyqWnTpm5XoEmTJq6zEAAktzFjxujBBx/UgQMHlDt3bjct2f7uAYCUdqnn5+xxhohhYRoUGOsb9NjKlTrep48LCo4ePRp7zHYCYjoJ3XzzzW7mAACkJEspsqnnbdq0cTVM1tq4e/fu6tu3ryIjI71eHgAkiB2DEHBK0lXWXUdhzKYPlyhhlcWuRiCmeLhatWq0DwTgCWt3/Oyzz+q9995z96tUqaJRo0apaNGiXi8NQIgqfYnn5/ReDAE/JDUosG47l/oP1Cuv+F5n+HB5qnhxPTJwoNasWaPNmzfrrbfecgXFBAUAvGI7lAMGDHCpjTlz5tTSpUtVrlw5V/MEAIGIwCAExG22Gb6qdumiMmXK0F4UQECxuqbVq1e7HYN9+/a53cxu3brFGaAIAIGAwCCcA4Off5Zmz1aoIEACEKgsfcgGKFqtgbGdBOtaZDucABAoCAzC+YT42mtdCk6oIDAAEMisTXL//v01adIk163ICpPLly+vL7/80uulAYBDYBDkrBORb5xOMtYYTJ0qNWgg5colZcwoXXONTeuR9u07/+v98IPUooXveVmySLVrS4sXn/tzVo9g7231Cb/9JrVvL112mZQpk1SxojRpUpJ+nTWSTibpmQCQeqxLkdVD2Y6BtTS1LkaPPvqom4MAAF4iMAhyVnOerFmqfftK1m973jypQgVLjrXhANIbb1hLDekvm5gQjxUrpKpVpa1bpUaNfB2CFiyQ6tWTfvwx/ufYz1aqJC1b5vu5cuWklSt97zlz5kUv/Yik9Rf9LABIfTb0bO7cuXrhhRfcfRuKZl3UNm7c6PXSAIQxAoMgl8Apd9IsXy699JKUNau0cKE0a5Y0apS0aZPUqpW0YYPUpUv8z/3vf31Bxfff+56zZo3UrZuNApX69Yv/OZ9+KnXo4Htde47tLrzzjnT6tNSnT5I7NAFAMLCZBq+99pqmT5+uvHnzul0ESy2yackA4AUCgyB3IDlf7P33fSfljz/u2x2IYUPB7Jil+owbJ23ffu5za9SQunaN+5gFGcZ2DuJTrJj02mtSmjP+M3zsMV8q0pIlUlTURf8KBy/6GQDgrUaNGmnt2rWqXbu2Dh06pPbt27vJyWcOawSA1EBgEOSSNSP12299X++669xj+fJJDRv6AodFi849bsfOliePlDu3tHNn/O9Xp45V48V9LF06X8Bw4oS0e/dF/wpk6AIIRpdffrlmzZqlnj17upbLQ4cOVeXKlfWzdY8DgFRCYBDkkvVEeMcO39eEhp7FPP7HH+ceK1Qo/udky5bwlf/zPcckocc319cABKt06dKpd+/emjlzpvLnz68ff/xRFStW1KeWdgkAqYDAIMil6h/g+QaHnZkOlFhJec4FMOcYQLCrX7++qzeoV6+ejhw5onvvvdfdDh8+7PXSAIQ4AoMglzE5X+zyy31ft21LuIuQueIKhcXnAQAeKVCggGbMmOF2ENKkSeN2DSpVquR2EQAgpRAYBLlMyfliN93k+xpfR4y//5ZmzPDtGlihcTh8HgDgobRp07qagzlz5rgaBKs3sOBg2LBhio6O9np5AEIQgUGQuyw5X8xakVp6z3vv+eYSxLAaAetUZB0ybr/dGnArUOX1egEAkMysW5GlFln3IhuC1qlTJ3Xo0EEHD9KHDUDyIjAIcuWS88UqV5ZefVU6cECqVs03/bhdO+nqq6XRo31Dy2xeQQAr7/UCACAFXHbZZZo6dar69u3rdhJGjhzpCpOtzSkAJBcCgyBnfYJyJWdBcY8e0uTJdonKN/Bs7FjfHINnn5WWLpXy51cg754k0OcIAIKe1Ro8//zzmj9/vgoVKqQNGzaoSpUq+vDDD0ktApAsIqKT+LdJ6dKl3dd169Ylz0qQZDZB4JuLfZKlBWXOLJUqZX+ICgVNJE31ehEAkAp2797tOhVNtgs5klq3bq2PPvpIOXLk8HppADx0qefn7BiEgApJeVJMDUHx4grrzwEAglCePHk0ceJE9e/f380/GDNmjCpUqKCVK1d6vTQAQYzAINxOiFetklq1kpo3T3jKcZAiMAAQTmxCcvfu3fXtt9+qSJEi2rx5s6pXr66BAweSWgQgSQgMQkCli/nh336Txo2TcueW+veX2rRRqKjo9QIAwANVq1bV6tWr1bJlS0VFRalr16664447tHfvXq+XBiDIEBiEgCKSqif2h1u2lE6elLZskbp3V6ioTeExgDCWK1cujR07VgMGDFBkZKTGjRun8uXLa6k1jQCARCIwCBGPKryF++8PAJZaZLsFixcv1lVXXaWtW7eqZs2aevvtt0ktApAoBAYh4s4wHu5VwDZCvF4EAAQIm2+watUqtWrVSidPnnR1CC1atHCdjADgfAgMQkQGSZ0Unh6UlN7rRQBAALG2paNHj9agQYOUIUMG19a0XLlybjcBABJCYBBCHratZIWXtJIe8noRABCgqUWdO3fWkiVLVKJECW3fvl21atXSG2+8odOnT3u9PAABiMAgxKYgN1N4aUHRMQCcV9myZd18g/bt2+vUqVNuenLTpk31999/e700AAGGwCDEvOFPKwoHGSW97vUiACAIZMuWTSNGjNDQoUOVMWNGTZ8+3QUMCxYs8HppAAIIgUGIuU5Sb4WH1ySV9HoRABBEqUUPPPCAli9frmuvvVY7duxQ3bp11adPH7eTAAAEBiHIphNUUWirIamr14sAgCB0/fXXa8WKFerYsaOrNejZs6caNWqkv/76y+ulAfAYgUGIFuR+EsIpRZZC9LH/9wQAXLwsWbJo+PDh7pY5c2bNnj1bZcqUcV8BhC8CgxAVyilFpBABQPKwXQNLLbJdBNsxaNCggXr16kVqERCmCAxCPKWonkJLI1KIACBZlSpVSkuXLlWnTp3chOTevXurfv36rgYBQHghMAhhlmozTlJlhYZqkr4mhQgAkp2lEw0ZMkQjR45U1qxZNW/ePNe1aMaMGV4vDUAqIjAIcdkkTbNiMwW3GyVNsbxYrxcCACHMZh3YzAOrN7A5B40bN1aPHj108uRJr5cGIBUQGISB3JLmSCqv4FTJv/5cXi8EAMJAyZIl3bRkm5ps+vbt69qa/v77714vDUAKIzAIE5f5T65rKbjUlWQ9MvJ4vRAACCM2BG3QoEEaPXq0smfProULF7rUoilTbO8WQKgiMAgjOSRNl/SYAl+EpCckTfWnQwEAUl/r1q21atUqVahQQbt371azZs30zDPP6MSJE14vDUAKIDAIM5kkDZQ0V1IxBabikuZLetc/swAA4J3ixYtr0aJF6trV1xOuf//+uummm7Rt2zavlwYgmREYhKk6kr4PsN2DmF0CW9dNXi8GABArQ4YMGjBggMaNG6ecOXO69qaWWjR+/HivlwYgGREYhLGsZ+we2FV6L5U4Y5cgs8drAQDEr2XLllq9erWqVKmiffv26bbbblO3bt0UFRXl9dIAJAMCA7jdg/X+mQcNPBhYNkHSz+wSAEBQKFq0qBYsWKDu3W2MptxOQo0aNbRlyxavlwbgEhEYwElnV4IkzZT0i6QnJeVMofeytqP2z8lGfzF0C4aWAUBQSZ8+vas1mDRpknLnzq0VK1aoXLly+uqrr7xeGoBLQGCAc5SU9LakPyR9KulBSeUkRSbx9SL9MxQekvQ//+v2l3R1Mq8bAJC6rEvRmjVr3I7BgQMH1KpVK3Xp0kXHjh3zemkAkiAiOjo6OilPLF26tPu6bt26pDwdQei4pB8krfTfbNP4qKRj/mMZ/F2EMvlrFir4b9f7jwEAQpO1L3355Zf1+uuvu/tWmDxmzBiVKGEVZABSy6WenxMYAACAZDF9+nR16NBB//zzj7JmzaqPPvpI7dq183pZQNgofYnn56QSAQCAZNG4cWOtXbtWtWvX1qFDh9S+fXs99NBDOnrU9pcBBDoCAwAAkGwuv/xyzZo1Sz179lRERISGDBmiypUra/16638HIJARGAAAgGSVLl069e7dWzNnzlT+/Pn1448/qkKFCvrss8+8XhqA8yAwAAAAKaJ+/fqua1G9evV05MgRdezYUffdd58OHz7s9dIAxIPAAAAApJgCBQpoxowZbgchTZo0Gj58uEstonkJEHgIDAAAQIpKmzatqzmYM2eOChYsqJ9++kmVKlXSsGHDlMTmiABSAIEBAABIFdatyFKLGjVq5DoVderUybU3PXjwoNdLA0BgAAAAUlO+fPk0depU9e3b1+0kjBw5UhUrVnRtTgF4i8AAAACkKqs1eP755zV//nwVKlRIGzZsUJUqVfThhx+SWgR4iMAAAAB4okaNGi61qFmzZjp+/Lg6d+6stm3b6sCBA14vDQhLBAYAAMAzefLk0cSJE9W/f383/2DMmDEqX768Vq5c6fXSgLBDYAAAADxlE5K7d++ub7/9VkWKFNHmzZtVvXp1DRw4kNQiIBURGAAAgIBQtWpVrV69Wi1btlRUVJS6du2qO++8U/v27fN6aUBYIDAAAAABI1euXBo7dqwGDBigyMhI9325cuW0bNkyr5cGhDwCAwAAEHCpRbZbsHjxYl111VXaunWrK1R+++23SS0CUhCBAQAACEg232DVqlVq1aqVTp486eoQbr31Vu3Zs8frpQEhicAAAAAErBw5cmj06NEaNGiQMmTIoEmTJqls2bJuNwFA8iIwAAAAAZ9aZDMOlixZohIlSmj79u2qVauW3njjDZ0+fdrr5QEhg8AAAAAEBdspsPkG7du316lTp9z0ZBuO9vfff3u9NCAkEBgAAICgkS1bNo0YMUJDhgxRxowZNW3aNBcwLFiwwOulAUGPwAAAAARdalGnTp1cC9Nrr71WO3bsUN26ddWnTx+3kwAgaQgMAABAULrhhhu0YsUKdezY0dUa9OzZU40bN9Zff/3l9dKAoERgAAAAglaWLFk0fPhwd8ucObNmzZqlMmXKaM6cOV4vDQg6BAYAACDo2a7B8uXLdf3117sdg/r166tXr16kFgEXgcAAAACEhFKlSmnp0qWu/sAmJPfu3dsFCFaDAODCCAwAAEDIsHQi61g0cuRIZc2aVfPmzXNdi2bOnOn10oCAR2AAAABCjs06sJkHVm9gcw4aNWqkHj166OTJk14vDQhYBAYAACAklSxZ0k1LtqnJpm/fvq6t6e+//+710oCARGAAAABClg1BGzRokEaPHq3s2bNr4cKFLrVoypQpXi8NCDgEBgAAIOS1bt1aq1atUoUKFbR79241a9ZMzzzzjE6cOOH10oCAQWAAAADCQvHixbVo0SJ17drV3e/fv79q1aqlbdu2eb00ICAQGAAAgLCRIUMGDRgwQGPHjlXOnDldDYKlFk2YMMHrpQGeIzAAAABh57bbbtPq1atVuXJl7du3Ty1btlS3bt0UFRXl9dIAzxAYAACAsFS0aFF9++236t69u7tvOwk1atTQli1bvF4a4AkCAwAAELbSp0/vag0mTpyo3Llza8WKFSpXrpy++uorr5cGpDoCAwAAEPaaN2+uNWvWuB2DAwcOqFWrVurSpYuOHTvm9dKAVENgAAAAIKlw4cKaO3eunn/+eXff5h9Uq1ZNGzdu9HppQKogMAAAAPCLjIx0E5KnTZumvHnzul2E8uXL64svvvB6aUCKIzAAAAA4S+PGjV1QYHMODh06pPbt2+uhhx7S0aNHvV4akGIIDAAAAOJxxRVXaPbs2erZs6ciIiI0ZMgQValSRevXr/d6aUCKIDAAAABIQLp06dS7d2/NnDlT+fPn1w8//KAKFSros88+83ppQLIjMAAAALiA+vXru9Sim2++WUeOHFHHjh1133336fDhw14vDUg2BAYAAACJUKBAAbdzYDsIadKk0fDhw93k5HXr1nm9NCBZEBgAAAAkUtq0aV3NwZw5c1SwYEH99NNPqlSpkoYNG6bo6GivlwdcEgIDAACAi1S7dm2XWtSoUSPXqahTp07q0KGDDh486PXSgCQjMAAAAEiCfPnyaerUqW7uge0kjBw5UhUrVtTatWu9XhqQJAQGAAAASWS1BjYped68eSpUqJA2bNjgWpoOHjyY1CIEHQIDAACAS1SzZk2XWtS0aVMdP35cjzzyiNq2basDBw54vTQg0QgMAAAAkkGePHk0ceJE9e/f380/GDNmjMqXL69Vq1Z5vTQgUQgMAAAAkjG1qHv37vr2229VpEgRbd68WdWqVdPAgQNJLULAIzAAAABIZlWrVtXq1avVsmVLRUVFqWvXrrrzzju1b98+r5cGJIjAAAAAIAXkypVLY8eO1YABAxQZGem+L1eunJYtW+b10oB4ERgAAACkkIiICLdbsHjxYl111VXaunWratSoobfffpvUIgQcAgMAAIAUZvMNrAjZ0olOnjzp6hBuvfVW7dmzx+ulAbEIDAAAAFJBjhw5XKeiQYMGKUOGDJo0aZLKli3rdhOAQEBgAAAAkIqpRZ07d9aSJUtUokQJbd++XbVq1VK/fv10+vRpr5eHMEdgAAAAkMpsp2DlypVq3769Tp06peeee07NmjXT33//7fXSEMYIDAAAADyQLVs2jRgxQkOGDFHGjBk1bdo0FzAsWLDA66UhTBEYAAAAeJha1KlTJ9fC9Nprr9WOHTtUt25d9enTx+0kAKmJwAAAAMBjN9xwg1asWKGOHTu6WoOePXuqcePG+uuvv7xeGsIIgQEAAEAAyJIli4YPH+5umTNn1qxZs1xq0Zw5c7xeGsIEgQEAAEAAsV2D5cuXq3Tp0vrzzz9Vv3599erVi9QipDgCAwAAgABTqlQpV3dg9Qc2Ibl3794uQLAaBCClEBgAAAAEIEsnso5FI0eOVNasWTVv3jyXWjRz5kyvl4YQRWAAAAAQwGzWgc08KFOmjJtz0KhRI/Xo0UMnT570emkIMQQGAAAAAa5kyZJuWrJNTTZ9+/Z1bU1///13r5eGEEJgAAAAEARsCNqgQYM0evRoNxxt4cKFLrVo6tSpXi8NIYLAAAAAIIi0bt1aq1evVoUKFbR79241bdpUzz77rE6cOOH10hDkCAwAAACCTPHixbVo0SI9/vjj7v6bb76pWrVqadu2bV4vDUGMwAAAACAIZciQQe+9957Gjh2rnDlzuhqEcuXKacKECV4vDUGKwAAAACCI3XbbbS61qHLlytq7d69atmypbt26KSoqyuulIcgQGAAAAAS5okWL6ttvv1X37t3d/QEDBqhGjRrasmWL10tDECEwAAAACAHp06dX//79NXHiROXOnVsrVqxwqUVfffWV10tDkCAwAAAACCHNmzfXmjVrVL16dR04cECtWrVSly5ddOzYMa+XhgBHYAAAABBiChcurHnz5un55593923+gQUKGzdu9HppCGAEBgAAACEoMjLSTUieNm2a8ubN6wqUy5cvr1GjRnm9NAQoAgMAAIAQ1rhxY5daZHMODh06pHbt2unhhx/W0aNHvV4aAgyBAQAAQIi74oorNHv2bPXs2VMRERH66KOPVKVKFa1fv97rpSGAEBgAAACEgXTp0ql3796aOXOm8ufPrx9++EEVKlTQZ5995vXSECAIDAAAAMJI/fr1XWrRzTffrCNHjqhjx4667777dPjwYa+XBo8RGAAAAISZAgUKuJ0D20FIkyaNhg8f7iYnr1u3zuulwUMEBgAAAGEobdq0rubAag8KFiyon376SZUqVdLHH3+s6Ohor5cHDxAYAAAAhLE6deq41KKGDRu6TkUPPPCAOnTo4DoYIbwQGAAAAIS5fPnyuXkHNvfAdhJGjhzpCpPXrl2rQGd7G8cl7ZO0S9J+SVH+x3FxCAwAAADgag1sUrJNTC5UqJA2bNjgWpoOHjw4YFKLdkv6RlJfSXdKKmYpUZIySsolKb+knJIyWBcmSddKukvS25LmSzrg9S8Q4CKik/gnXbp0afeVIhUAAIDQsnv3btetaMqUKe5+69atNWTIEGXPnj1V1xHtP6EfKmmRpK2X+HoRkkpY+pSkhyWVV2gpfYnn5+wYAAAAII48efJo4sSJ6t+/v5t/MGbMGJUvX16rVq1Klfe3dKD37URXUl1JI5MhKIgJNDZI+khSBUlVJdkUh2PJ8NqhgMAAAAAA8aYWde/eXd9++62KFCmizZs3q1q1anr//ffjTS0aPXp07A5DUv0o6RGb1CzpcUk/K2UtldTR/37PJFPwEcwIDAAAAJCgqlWravXq1WrZsqWioqL0+OOP684779S+fVbu+29Q0LZtW7Vo0UIbN2686PewK/bPSSojabCk1B61tkdSf0nXSXpL0imFJwIDAAAAnFeuXLk0duxYDRgwQJGRke77cuXKadmyZdq0aZMefPBB93OnT5/Wq6++etFX7ctJ6mfPl7eOSXpa0k2SflH4ofgYAAAAibZixQq1adNGW7ZscfUH+fPn1x9//BEnBcmGpV1zzTUXPAnv5b9S73VAEJ+MkvpI6ubvfBQMKD4GAABAqqlYsaIrQrZ0opMnT8YJChK7a7AugHYJErN78LvCA4EBAAAALkqOHDlcTUFCvvjiC61fvz7B1CE72Y7/aOD5TlJNSZsU+ggMAAAAcFF+/fVXPfDAAwket12D3r17n/P4HEn1JO1VcNnmDw6+V2gjMAAAAMBF6dChg/bvt2kDCRs1apSrNYixUFIzDzoOJZe/JNUP8aJkAgMAAAAk2qlTp1wnogux/jbW2tSsltRU0lEFt78lNZD0m0JTOq8XAAAAgOCRNm1aff/995o/f75++OEH973dLL3obHPmzNGCDRvUqmRJHVBo2C6pkaSVkjIrtNCuFAAAAJfs4MGD+vHHH12QYF2LJk2apFOnT6vi9u2aGhmpUNNN0jsKLJd6fk5gAAAAgBTxuaS7FJoiJC3wFyUHCuYYAAAAIOD8KclXYRCaoiXdJ+mIQgeBAQAAAJL9pPkRSXsU2jZJelGhg8AAAAAAyWqUpAkKDwP8rVhDAYEBAAAAks1pSb0UXrsj/1FoIDAAAABAspktaaPCy6wQGXxGYAAAAIBkM0jh6UMFPwIDAAAAJNvwr4kKT59IOqzgRmAAAACAZDHYX2MQjvb7i66DGYEBAAAALlmUpCEKb4MU3AgMAAAAcEFbt25VRESE6tSpE+/xZZJ2KUBt3SpFREgJrD25rJL0u4IXgQEAAAAu2UqvFxAgVip4ERgAAAAgrE+Ik9NKBS8CAwAAAIT1CXFyWqngRWAAAACAi3LgwAE98cQTKly4sDJmzKhrrrtOP73zjnT6rJ5EltdftGj8LzJ8uO/4K6/EfdzqAOxxqwv4/HOpalUpWzYpZ85/f+bwYen5532vnTGjdPXV0quvSidO+B6z5yfGsWPSsGHSrbdKV10lZcrke59ataRRCfQYuvde3+vPmyfNmCHVret7jj22c6emRka6z+XUqVPxPv3zzz93tRodO3ZUoCEwAAAAQKIdP35cN998sz777DNVrlxZDRo00LZt26SnnpLuvz/53qhvX6lDByl9eqlZM+n662MWINWvL73xhnTwoO/YdddJr78utW59ce9hwUenTtKKFb6AwgKEsmWlJUukdu3ODVrOZEFLkya+IMW+VqokZc4stWih33//XdOnT4/3aUOG+Ho3PfTQQwo06bxeAAAAAILHkiVLdOONN2rjxo3Kmzeve6z35s3qZVfZP/1UatnSd7tUn30mzZkj1a4d9/E33/SduFeu7LtiH7OTYCf5tobtNmYtkS67TPrmG6levbi7DL/+Kt18s28XwnYI4tv1sBN821Vo0ybu4488Io0d6wKApk2bxjm0adMmzZs3T9ddd51q1KihQMOOAQAAAC5K//79Y4MCk6F4calnT9+d999Pnjd54IFzgwLz4Ye+r2+9FTe9yE7eX3754t4jTx7f7sPZqUfFikkvvuhLjZo0Kf7n2kn/2UGBqV9fBa6+WlOmTNHOnTvjHBo6dKj7+uCDDyoQERgAAAAg0XLnzu3Sh8501P7HUm/M4sXn1hokRYsW5z5mKUt//CEVKCDVrHnu8fhO1BNj4UKpTx+pc2fpvvt8uwRffuk7tnFj4tdnIiJU66GHdPLkSX3yySexD584cULDhw9XhgwZdM899ygQkUoEAACARCtSpMg5jx2z/8mRw3cFf98+ae9e39X4S3Hllec+FnMFvnDh+J8TU6Rsa0iM/ful22/3pSwlxOoYErs+v4r33acJPXtq2LBheuGFF1yx8aRJk/TXX3+pXbt2ynOpn00KYccAAAAAlySRPYDiutCugnUbSmnPPfdvHcO8edI//0gnT0rR0b76BWPfX+T6sufNqzvuuENbtmzR7NmzgyKNyBAYAAAAINF+++23cx7LZP9z4IDvSn1My08TGSkdOhT/C11MkXCMggXP/1y7up/Y3QIzbpyUNq00caIvOLAr+XbfbNmipLLP4xErQvZ3IbLPbMaMGSpRooTqWnvTAEVgAAAAgETbvXt37FXwGO7aeUzf/2rV/j25thP53bt9t7PNmnXxb25pTFdcIf35p6+W4WwxdQGJZSlP2bP7bmcbM0ZJZZ/HTTfdpNKlS2v8+PHq16+fTp8+rU7WGjWAERgAAADgojz99NMuQIgRbe09e/f23enS5d8fjOkqZIW9Z+rXz1fwmxT+K/Hq3t1XI3BmYXLMGhKrZElfcDB6dNzHbVjb3LlJW5+kmH5NDz/8sKKiovTf//5XkZGRuteKmgMYgQEAAAASrWrVqkqTJo2uvvpql0ffokULvWLDx6xb0N13+4p5z8zht9Sid9+VypWT7rxTuuYa3+CwRx9N2gKeecY3DdlmGVib1FatfB2CSpeWypTxFQVbClNivPCC72vbtr4ZCO3b+17n6aelJ59M2voklfV/te5DmW3omWx22q3Kly+fAhmBAQAAABLN2m3OmTNH7du3d8POLHf+ysKFlbZ/f2n48Lg/bCfZVtxbp460YYNvmJidzH/3nW9ScNIW4HudZ5+VsmTx1QesW+fbQbAr/3/9lfiOSHfdJU2Z4gs01qyRpk2TLr/ct+aE2pFewFXW0tX/fY4cOVS+fPmALzqOEREdnVCp9flZzpRZZ38QAAAACGuVJS33ehG2i2A1Do0b+07yPdDKyhP832/fvl3FihVT4cKFXYcia1uaki71/JwdAwAAAFyyCqn5ZqtXn9vu1LoIPfyw73tLaQqAz+H111/XqVOn1KVLlxQPCpIDA84AAAAQXIGBTVm2wuMbbvClDVkL1ZUrpePHfSlAVivgkXy//KJOb76pX3/91aVcFSpUKLZ1aaAjMAAAAMAlS2LFQNI8/rivnej330t79viGjZUt69spsJNwj67OR0jKs3Onm3icKVMm1apVSwMHDlTWrFkVDKgxAAAAwCWzE0o7O/xZ4auJpKkevj81BgAAAPCcXS1PYgPSkPGoghuBAQAAAJJFB0lZFJ6K+HcMghmBAQAAAJJFDmsIpPD0iKS0Cm4EBgAAAEg2nRV+0ku6X8GPwAAAAADJpoykGgovraxNqYIfgQEAAACSVf8wOsnMIulVhYZw+TMDAABAKqkqqbvCw5uSiik0EBgAAAAg2f1H0jUKbXUlPazQQWAAAACAZJdJ0vAQPtnMImlYiP1+ofS7AAAAIICEckrRmyGUQhSDwAAAAAApmlJUXqGlWYilEMUgMAAAAECKphRNl3StQsNNkkaH6El0KP5OAAAACCCXSfpGUhEFt3KSJknKrNBEYAAAAIAUV0jSfElXKzhVkTRLUg6FLgIDAAAApArbMVgo6UYFl3r+oCC3QhuBAQAAAFJNfv/OQTsFvghJT0iaLCmrQh+BAQAAAFJVTkmfSxorKZ8CU3F/APOupIwKDwQGAAAA8MRtkn4KsN2DCP8uwff+DkThhMAAAAAAnslzxu6BpRl56eozdglCtfPQ+RAYAAAAICB2D36VNFxSZQ+Ki7+W9HMY7hKcicAAAAAAATMMraOkpZKWS7o/BfP7re1oN0nr/R2HbpeUTuEt3H9/AAAABKCKkoZJetOfarRI0kpJGy/havh1kipIqiOptaQsybzmYEdgAAAAgIBlswMe89/Mfkmr/EGC3bZKOibpqKQoSRn8Ow+201DCH2BYMFCGQOCCCAwAAAAQNCwFqK7/huRFjQEAAAAAAgMAAAAABAYAAAAACAwAAAAAGAIDAAAAAAQGAAAAAAgMAAAAABAYAAAAADAEBgAAAAAIDAAAAAAQGAAAAAAgMAAAAABgIqKjo6OT8lFky5ZNJ06cUPHixfkkAQAAAI9t3rxZkZGROnjwYOruGGTJksW9MQAAAADv2bm5naOn+o4BAAAAgNBBjQEAAAAAAgMAAAAABAYAAAAACAwAAAAAGAIDAAAAAAQGAAAAAAgMAAAAABAYAAAAADAEBgAAAAAIDAAAAAAQGAAAAAAgMAAAAABgCAwAAAAA4f8Br2jXDEmmscMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 960x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# ===== 代码块3：定义贝叶斯网络相关类 =====\n",
    "class BayesNode(Node):\n",
    "    def __init__(self, name):\n",
    "        super().__init__(name)\n",
    "        self.cpt = {}\n",
    "    \n",
    "    def cpt_score(self, evidence):  # 只返回为1的概率\n",
    "        # 我们假设，调用这个方法的时候，该节点的所有支撑变量都进了evidence(这需要推理过程确保)\n",
    "        key_tuple = []\n",
    "        for parent in self.parents:\n",
    "            if parent in evidence:\n",
    "                key_tuple.append(evidence[parent])\n",
    "        key_tuple = tuple(key_tuple)\n",
    "        # 对于先验概率，其tuple应该是一个空值，这样可以直接算出需要的p\n",
    "        p = self.cpt[key_tuple]\n",
    "        return p\n",
    "    \n",
    "    def __repr__(self):\n",
    "        return \"{}{}\".format(self.name, tuple(p for p in self.parents))\n",
    "\n",
    "class BayesNet(Graph):\n",
    "    def __init__(self, name='', cpt_pd=\"./burglary_cpts.csv\"):\n",
    "        super().__init__(name)\n",
    "        self.cpt_pd = pd.read_csv(cpt_pd)\n",
    "        self.parse_cpts(self.cpt_pd)\n",
    "    \n",
    "    def parse_cpts(self, data_pd):\n",
    "        # 首先加载各个点\n",
    "        # 注意这里必须严格确保源文件中列的顺序，因为检查的时候是从先验概率一层一层补充检查的\n",
    "        for col in data_pd.columns[2:]:\n",
    "            self.add_node(BayesNode(col))\n",
    "        \n",
    "        # 构造各个点内部关系\n",
    "        for idx, row in data_pd.iterrows():\n",
    "            p, key = row['p'], row['key']\n",
    "            row['p'], row['key'] = 0, 0\n",
    "            row = row.where(row != 0).dropna()\n",
    "            row_node = self.node_map[key]\n",
    "            for col, val in row.items():\n",
    "                row_node.parents[col] = 1\n",
    "            vals = tuple(val for _, val in row.items())\n",
    "            row_node.cpt[vals] = p\n",
    "    \n",
    "    def ask(self, x, evidence, print_formula=False):\n",
    "        result = {1: 0, -1: 0}\n",
    "        if x in evidence:  # 如果x的状态已知\n",
    "            result[evidence[x]] = 1.0\n",
    "        else:  # 拿出来所有的node，逐一排查\n",
    "            cal_seq = []\n",
    "            cal_seq.append(\"1:\")\n",
    "            result[1] = self.recurse_ask(self.node_list, {**evidence, x: 1}, cal_seq)\n",
    "            cal_seq.append(\"\\n-1:\")\n",
    "            result[-1] = self.recurse_ask(self.node_list, {**evidence, x: -1}, cal_seq)\n",
    "            result_sum = result[1] + result[-1]\n",
    "            result[1] = result[1] / result_sum\n",
    "            result[-1] = result[-1] / result_sum\n",
    "        \n",
    "        if print_formula:\n",
    "            print(\"\".join(cal_seq))\n",
    "        return result\n",
    "    \n",
    "    def recurse_ask(self, check_list, evidence, cal_seq):\n",
    "        if not check_list:\n",
    "            cal_seq.append(\"1\")\n",
    "            return 1.0\n",
    "        \n",
    "        check_node, rest = check_list[0], check_list[1:]\n",
    "        if check_node.name in evidence:  # 如果这个状态是一个已知的值\n",
    "            cpt_score = check_node.cpt_score(evidence)\n",
    "            if evidence[check_node.name] == -1:\n",
    "                cpt_score = 1 - cpt_score\n",
    "            cal_seq.append(\"({}*\".format(cpt_score))\n",
    "            recurse_val = self.recurse_ask(rest, evidence, cal_seq)\n",
    "            cal_seq.append(\")\")\n",
    "            p = cpt_score * recurse_val\n",
    "        else:  # 如果这个状态未知，就需要计算+-1的所有情况\n",
    "            cpt_score = check_node.cpt_score(evidence)\n",
    "            cal_seq.append(\"({}*(\".format(cpt_score))\n",
    "            p = cpt_score * self.recurse_ask(rest, {**evidence, check_node.name: 1}, cal_seq)\n",
    "            cal_seq.append(\")+{}*\".format(1 - cpt_score))\n",
    "            p += (1 - cpt_score) * self.recurse_ask(rest, {**evidence, check_node.name: -1}, cal_seq)\n",
    "            cal_seq.append(\")\")\n",
    "        return p\n",
    "\n",
    "burglary = BayesNet('burglary')\n",
    "burglary.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "845fe7c9-cbeb-4dfa-b0b9-f60aaf7e2538",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "burglary {1: 0.001, -1: 0.999}\n",
      "earthquake {1: 0.002, -1: 0.998}\n",
      "alarm {1: 0.0025164420000000002, -1: 0.997483558}\n",
      "john {1: 0.0521389757, -1: 0.9478610243}\n",
      "mary {1: 0.011736344980000002, -1: 0.98826365502}\n"
     ]
    }
   ],
   "source": [
    "# ===== 代码块4：创建贝叶斯网络并进行查询 =====\n",
    "burglary = BayesNet('burglary')\n",
    "\n",
    "# 计算贝叶斯网络中的独立概率\n",
    "print('burglary', burglary.ask('burglary', {}))\n",
    "print('earthquake', burglary.ask('earthquake', {}))\n",
    "print('alarm', burglary.ask('alarm', {}))\n",
    "print('john', burglary.ask('john', {}))\n",
    "print('mary', burglary.ask('mary', {}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a880802e-def8-4208-a082-cb4bfbcfd4da",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1:(0.001*(0.002*((0.95*((0.9*((0.7*(1)+0.30000000000000004*1))+0.09999999999999998*(0.7*(1)+0.30000000000000004*1)))+0.050000000000000044*(0.05*((0.01*(1)+0.99*1))+0.95*(0.01*(1)+0.99*1))))+0.998*(0.94*((0.9*((0.7*(1)+0.30000000000000004*1))+0.09999999999999998*(0.7*(1)+0.30000000000000004*1)))+0.06000000000000005*(0.05*((0.01*(1)+0.99*1))+0.95*(0.01*(1)+0.99*1)))))\n",
      "-1:(0.999*(0.002*((0.29*((0.9*((0.7*(1)+0.30000000000000004*1))+0.09999999999999998*(0.7*(1)+0.30000000000000004*1)))+0.71*(0.05*((0.01*(1)+0.99*1))+0.95*(0.01*(1)+0.99*1))))+0.998*(0.001*((0.9*((0.7*(1)+0.30000000000000004*1))+0.09999999999999998*(0.7*(1)+0.30000000000000004*1)))+0.999*(0.05*((0.01*(1)+0.99*1))+0.95*(0.01*(1)+0.99*1)))))\n",
      "{1: 0.001, -1: 0.999}\n"
     ]
    }
   ],
   "source": [
    "# ===== 代码块5：检查计算公式 =====\n",
    "# 如果要检查计算公式，注意最终的值取了归一化\n",
    "print(burglary.ask('burglary', {}, True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9ec1a18d-98ad-44da-9ed0-34c2825d279a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{1: 0.2841718353643929, -1: 0.7158281646356071}\n",
      "{1: 0.005129858133401302, -1: 0.9948701418665987}\n",
      "{1: 0.7521522838357262, -1: 0.24784771616427376}\n",
      "{1: 0.8999999999999999, -1: 0.09999999999999998}\n",
      "{1: 0.9, -1: 0.09999999999999999}\n"
     ]
    }
   ],
   "source": [
    "# ===== 代码块6：不同证据下的查询 =====\n",
    "# john和mary都打电话\n",
    "print(burglary.ask('burglary', {'john': 1, 'mary': 1}))\n",
    "\n",
    "# john打电话时，mary没有打电话\n",
    "print(burglary.ask('burglary', {'john': 1, 'mary': -1}))\n",
    "\n",
    "# 注意，同样的状态，但条件不一样的时候，概率也不相同\n",
    "# 下例和上例都有一个状态，即bljlm-1，但贝叶斯里条件和结论位置不同，不同状态的概率是不同的\n",
    "print(burglary.ask('john', {'burglary': 1, 'mary': -1}))\n",
    "\n",
    "# 对于john和mary，如果父类固定了，那么两者相互独立\n",
    "print(burglary.ask('john', {'alarm': 1}))\n",
    "print(burglary.ask('john', {'alarm': 1, 'mary': -1}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11fde9ee-d779-4050-9ef6-9bb1d53c689d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.21"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
